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N the fall of 1959, the Space Science Board of the 
National Academy of Sciences was informed of a 
proposal by the Lincoln Laboratory of the Massa- 
chusetts Institute of Technology thatjone or more belts 
of thin microwave dipoles be placed into orbit about 
the earth at an altitude of a few thousand kilometers 
as a means for achieving a reliable system of long-range 
radio communication. The Board was asked both by 
the Director, Lincoln Laboratory, and the government 
to consider the effects of this new communications 
technique from the standpoint of science as a whole 
and particularly to ascertain whether the first experi- 
mental test proposed by the Lincoln Laboratory would 
destroy or put any important limitation on activity 
in any field of basic scientific research. 

In response to this request, the Board appointed an 
ad hoc committee of scientists under,the chairmanship 
of Dr. O. G. Villard, Jr., to examine the consequences 
of the proposed experiment. Astronomers on this 
committee included G. M. Clemence, U. S. Naval 
Observatory ; C. H. Mayer, Naval Research Laboratory, 
and F. L. Whipple, Smithsonian Astrophysical Obser- 
vatory. The committee held several meetings, the last 
of which took place at the Lincoln Laboratory in 
April 1960 and was attended also by about a dozen 
optical and radio astronomers and by representatives 
of the Lincoln Laboratory and of several government 
agencies. 

After receiving the report of its ad hoc committee 
in June 1960, the Space Science Board arrived at a 
number of conclusions and recommendations: (1) The 
first exploratory test proposed by the Lincoln Labora- 
tory would probably not have an adverse effect on any 
branch of science. (2) The Board was seriously con- 

cerned with the harm which an operational system or 
systems might entail for optical and radio astronomical 


observations and strongly recommended that any 
planning for such a system or for large-scale future 
tests protect the interests of basic astronomical research 
and of science in general. (3) Full information on the 
scientific and operational aspects of the initial experi- 
ment should be published as soon as possible. (4) In 
view of the possible seriousness of the interference to 
radio astronomy, which is likely to result not only from 
dipole belts, but also from active communication 
satellites, protected frequency bands for radio 
astronomy should be established on a world-wide basis. 
Finally, the Board established a committee of radio 
astronomers to work closely with the developers of 
Project West Ford and to serve in an advisory capacity 
in connection with the question of interference to 
astronomy. 

The present composition of the liaison committee is 
as follows: John Findlay (National Radio Astronomy 
Observatory), F. T. Haddock (University of Michigan), 
William Liller (Harvard University), and A. E. Lilley 
(Harvard University). Dr. W. A. Hiltner (Yerkes 
Observatory) and Dr. Allan Sandage (Mount Wilson 
and Palomar Observatories) have recently been invited 
to join this committee. The primary task assigned to 
the committee was to review and to study carefully the 
theoretical predictions of the effects of the West Ford 
experiment on astronomy, which had also been under 
continual study by the Lincoln Laboratory, and to 
provide its scientific findings to the Space Science 
Board. Shortly after the committee began its work, the 
details of the West Ford experiment were described 
by Dr. W. E. Morrow, of the Lincoln Laboratory, in a 
paper presented before the URSI General Assembly in 
London in September 1960. In January 1961, the liaison 
committee met at the Lincoln Laboratory with Dr. 
Morrow and his associates. At this time it received an 
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up-to-date and thorough review of the status of the 
project and clarified a number of questions concerning 
the parameters of the experiment. 

The Space Science Board believes strongly that a 
discussion of proposed experiments that may have 
serious consequences for science generally should not be 
limited to committee consideration, however competent 
its members may be, but should be brought out for 
critical analysis by the scientific community as a whole. 
It has, therefore, recommended that the theoretical 
and technical aspects of the West Ford program be 
published in the astronomical literature together with 
a full report of the. studies of the laison group. The 
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three articles which follow are published in response to 
this recommendation. 

The Space Science Board wishes to express its warm 
appreciation for the readiness with which Dr. C. F. J. © 
Overhage, Director, Lincoln Laboratory and Dr. W. E. © 
Morrow have cooperated in furnishing the Board with 
all pertinent information regarding Project West Ford © 
and especially for their prompt agreement and assistance ~ 
in providing the details of the experiment for publi- — 
cation. The Board is also indebted to Dr. Dirk Brouwer, — 
Editor of the Astronomical Journal, for prompt publi- — 
cation of the papers in response to Space Science Board — 
arrangements with the American Astronomical Society. 
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Thin microwave dipoles orbiting in belts about the earth have been proposed for reliable long-range radio 
communication. This paper examines the possible effects of such microwave dipole belts on radio communi- 
cations, radar, and measurements by optical and radio astronomers. In addition, a comparison is made of 
the probabilities of micrometeorite and dipole collisions with space craft. Calculations have been carried out 
for a proposed experimental belt consisting of 35 kg of 8000 Mc dipoles orbiting at a few thousand kilometers 
altitude. It is shown that no interference would be caused to radio or radar services, that such a belt produces 
only a few percent change in the optical brightness of the night sky in the direction of the belt, and that the 
influence of the belt on radio astronomy observations is not measurable except when viewing the belt where 
it is illuminated by very high power radio equipment operating at frequencies near the dipole resonance. 
Investigation also shows that micrometeorite collisions with a space craft located in a dipole belt are an order 
of magnitude more likely than those caused by the dipoles. Furthermore, the small mass of the dipoles limits 


the collision effects to little more than surface scratches. 


INTRODUCTION 


OR the past 2} years a new method of providing 

reliable intercontinental communication has been 
under investigation at Lincoln Laboratory. This tech- 
nique involves the use of belts of orbiting microwave 
dipoles as an artificial microwave scattering medium. 
An initial study of this technique indicates that 
communication capacities of many tens, or perhaps 
hundreds of thousands, of bits per second can be 
achieved over distances of many thousands of miles 
(Morrow 1960). 

Figure 1 shows how two such belts might be used with 
highly directive microwave antennas to provide a 
multiplicity of high-capacity circuits having global 
coverage. The significant advantages of this technique 
over other satellite communication techniques are: 


(A) The belts are quite invulnerable to physical 
damage. 

(B) No complex electronic equipment is required in 
space. 

(C) Very large numbers of circuits can be provided 
by a single belt. : 

(D) Ground antennas are not required to track at 
high rates over large angles. 


Part of this study has been devoted to the possible 
problems that such a technique might produce because 
of the modification of the electromagnetic and physical 
properties of the space through which the dipoles orbit. 
This report is a summary of the findings of this study. 


PHYSICAL CHARACTERISTICS OF TEST AND 
OPERATIONAL BELTS 


Studies have been carried out to determine optimum 
dipole design from the point of view of maximizing 
the scattering cross section in the antenna beams for a 


~ given weight in orbit while at the same time minimizing 


* Operated with support from the U. S. Air Force. 


collision effects with space vehicles and electromagnetic 
effects at wavelengths other than that used for orbital 
scatter communication systems. While not complete 
in all respects, these studies indicate the desirability 
of dipoles resonant at wavelengths of 3 to 5 cm and 
having maximum geometric area-to-mass ratios of 
about 50 square centimeters per gram. The choice of a 
short wavelength permits the dipole mass to be reduced, 
thus minimizing collision hazards. Atmospheric and 
precipitation absorption effects prevent the use of still 
shorter wavelengths. The area-to-mass ratio is chosen 
sufficiently low to minimize perturbation of the belts by 
solar radiation pressure effects and micrometeorite 
collisions with the dipoles. The realization of a satis- 
factory area-to-mass ratio, while at the same time 
minimizing the dipole weight, leads to the use of high- 
density metals such as copper, tin, and tungsten which 
can be drawn into the thin diameters required. 


Fic. 1. Polar and equatorial dipole belts. 
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Fic. 2. Typical dipoles. 


Typical dipole designs using these metals are given 
in Table I. See Fig. 2 for typical copper dipoles 0.0025 
cm in diameter. 


It is planned to launch a test belt, using a piggy-back_ 


arrangement on an already programmed satellite. A 
payload of 100 lb on the vehicle is available. With this 
payload approximately 75 lb (35 kg) of dipoles can be 
launched into a belt. The rest of the payload is required 
for the dipole dispensing mechanism. A high-latitude 
inclination near-circular orbit with an altitude of a few 
thousand kilometers is planned. The dipoles will be 
released from a rotating dispenser with velocities 
ranging from about 0.8 to 3.0 m/sec. This will result 
in the cloud of dipoles spreading along the orbit ahead 
of and behind the dispenser. The cloud length will grow 
at the rate of about 2000 km per day. After about 30 
days the dipole belt will extend around the earth. The 
transverse dimensions of the belt are determined by 
differential perturbation of the dipole orbits. The first 
differential perturbation is that produced by the 
dispensing velocity distribution which will cause the 
belt to disperse about 32 km in a radial direction from 
the earth and about 8 km in a direction orthogonal to 
the orbit and a radial from the earth’s center. It is 
expected that other differential perturbations will occur 
because of the effects of solar radiation pressure, 
micrometeorites, and perhaps interaction between 
charged dipoles and the earth’s magnetic field. The 
extent to which these effects will disperse the belt are 
still under study. However, it is believed that the 
greater part of the dispersion will be less than about 
40 km in a radial direction at the end of 60 days and 
that probably it will be of the order of 100 km at the 
end of a year. (These figures are for an area-to-mass 
ratio of 50 cm?/g. The use of lower area-to-mass ratios 
will reduce the dispersion roughly in proportion.) 

An operational belt would probably be located at a 
higher altitude, perhaps 5000 to 6000 km, in order to 
obtain complete communication coverage and to 
minimize solar radiation pressure orbit perturbations 
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which are especially strong for polar orbits at altitudes 
of 3700 km due to a resonance phenomenon reported 
elsewhere (Jones and Shapiro 1960). Much larger 


payload weights could be used than in the case of the ~ 
test belt. The payloads will depend, of course, upon the ~ 


vehicles available, but it would appear that payloads 
of 1000 kg should be possible in the very near future and 
up to 3000 kg within a few years. 


An important aspect of the belt characteristics is — 


that of lifetime. At altitudes of greater than about — 
1000 km atmospheric friction cannot be depended upon ~ 
to bring the dipoles down. Depending upon the particu- © 


lar values of altitude and inclination that are realized, ” 
solar radiation pressure_effects can act to produce an_ 


ever increasing eccentricity which would bring perigee 
into the atmosphere within a period of one to two years. 


Alternatively, under other initial conditions of belt — 


altitude and inclination, the belt lifetime could be as 
much as one to two decades. Of course, due to the 
increasing dispersion of the belt, the dipole density per 


unit solid angle will have reduced enormously within — 


the first few years. In the case of operational belts at 
higher altitudes, the solar radiation pressure produces 
only small variations in the eccentricity and thus 
cannot at all be depended upon for limiting the lifetime 
of the dipoles in orbit. Since it is very desirable to limit 
the life of the dipoles to a few years for several reasons, 
development has been initiated on finite lifetime 
dipoles. A technique is being investigated which will 
result in the dipoles breaking up into small particles 
after a few years. This involves the transformation of 
white tin-alloy dipoles into gray tin powder under the 
temperature environment experienced in space. Solar 
radiation pressure will then bring the powder down into 
the earth’s atmosphere. Experiments have indicated 
that a successful finite lifetime dipole can be developed 
using this technique. 

Typical characteristics of the test belt are summarized 
below in Table II. 


ELECTROMAGNETIC PROPERTIES OF ORBITING 
DIPOLE BELTS 


The electromagnetic’ behavior of a dipole belt 
depends upon the density of the dipoles in the belt, the 
belt dimensions, and the electromagnetic scattering 
properties of the individual dipoles. The belt dimensions, 
together with typical dipole densities, are given in the 


TasLE I. Dipole designs. 


Dipole length 1.77 cm 0.7 inch 
Dipole diameter 
(For area-to-mass 
50 cm?/g) 
Tin 0.00350 cm 0.00138 inch 
Copper 0.00286 cm 0.00112 inch 
Tungsten 0.00131 cm 0.00052 inch 
Dipole mass 
Tin 0.00012 ¢g 
Copper 0.00010 g 
Tungsten 0.000047 g 
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Taste II. Test belt characteristics. 


Dipole resonance 8000 Mc 
Dipole length 1.77 cm 
Area-to-mass 50 cm?/g 
Belt altitude A few thousand kilometers 
Belt dimensions 8X40 km 60 days from launch 
Mass of dipoles in belt 35 kg 
Number of dipoles in belt 3.5108 
Average separation of 
dipoles 370 m 


Density of dipoles 21 per km? 60 days from launch 


previous section. The electromagnetic properties of a 
single dipole are discussed below. 

The dipoles consist of thin wires one-half wavelength 
long at the operating frequency with a diameter about 
1/1000 of the length. Such dipoles exhibit a directivity 
in their scattering pattern. The maximum response 
occurs with backscatter at resonance in a direction at 
right angles to the dipole length. A scattering value of 
0.89)? is obtained for this case This compares with 
nearly zero for backscattering in directions off the end 
of the dipole. In practice the dipoles will be tumbling 
randomly and thus on the average the scattering area 
will be somewhat less than 0.89\? per dipole. Calcu- 
lations have been carried out which show that under 
conditions of optimum polarization of the antennas 
0.159\? (2.22 sq cm at 8000 Mc) of scattering cross 
section is obtained on the average. Values as low as 
0.05\? (0.7 sq cm at 8000 Mc) result when the poorest 
choice of antenna polarization is made. 

The response at frequencies other than at the half- 
wave resonance of the dipole can be calculated from a 
knowledge of the impedance characteristics of thin 
dipoles. The results of such a calculation are shown in 
Fig. 3. It will be noticed that the maximum response 
occurs at resonance as might be expected. The cross 
section at lower frequencies falls off rapidly. Well below 
the resonance it decreases as the fourth power of the 
frequency (Rayleigh scattering region). The amplitude 
of the higher resonances falls off with frequency until 
the wavelength is comparable with the diameter of the 
dipole. This occurs in the near-infrared region. The 
cross section at shorter wavelengths corresponds to the 
projected geometric area of the dipole corrected by the 
albedo (reflectivity) of the dipole. The projected 
geometric cross section has a value of about 5x10-% 
sq cm and the albedo can range from 0.1 to 1.0 at optical 
wavelength depending upon surface conditions. A 
typical albedo value might be 0.3. 

From the above information typical dipole scattering 
cross sections can be deduced for various wavelengths. 
The scattering properties per unit volume in the test 
belt are found by using the above data, together with 
the information contained in the previous section on 
the density of dipoles in the belt. Table ITI summarizes 
the cross section versus wavelength of dipoles resonant 
at 4 cm wavelength. The cross sections are calculated 
‘on the basis of a dipole resonance cross section of 
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Fic. 3. Relative dipole scattering cross section. 


0.1? to allow for antenna polarizations other than 
optimum. 

The electromagnetic effects of the belt can now be 
determined. The areas of chief interest are: (1) the 
effects in the region of visible light, particularly as they 
might effect optical astronomical measurements; (2) 
the effects on radio astronomy measurements, both at 
the dipole resonance and at other frequencies; (3) the 
possibility of causing interference with communication, 
navigation and radar systems. 


OPTICAL WAVELENGTH EFFECTS 


The first effect to be considered is that of the attenua- 
tion experienced by visible light in traversing the belt. 
Examination of Table III shows that this will be 
extremely small for the case of the test belt: The 
maximum interception of light will be approximately 
the product of the belt thickness and the dipole density. 


40 kmX3X10~ cm?/km'= 1.20 cm?/km?. 


This is about 1 part in 10”. Clearly the light intensity 
from the sun, moon, planets, and stars will not be 


Taste III. Cross section vs wavelength of dipoles 
resonant at 4 cm wavelength. 


Cross section Cross section 


per dipole per cubic 
(Random polari- kilometer test 
zation of belt (60 days 
Wavelength Frequency antennas) after launch) 
4000 m 75 ke 6.31077! cm? 1.210 cm? 
400 m 0.75 Mc 6.31077 cm? 1.210715 cm? 
40 m 7.5Mc 6.31078 cm? 1.2107" em? 
4m 75 Mc 6.3107 cm? 1.21077 cm? 
40 cm 750 Mc 6.31075 cm? 1.2X10-? cm? 
10 cm 3300 Mc 4.0107 cm? 8) 5<1074) cm 
4cm 7500 Mc 1.6 cm? 32 cm? 
1 cm 33000 Mc 8 X107? cm? 1.6 cm? 
4mm 75000 Mc 2 X107 cm? 0.4 cm? 
Visible light and 
near infrared 1.510 cm? 30 <10-F em 
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affected to a measurable degree nor will there be any 
measurable angular scintillation effects. 

The second effect to be considered concerns the 
effect of scattered sunlight on nighttime observation of 
weak astronomical light sources, such as distant stars 
and galaxies, air glow, and zodiacal light. For the case 
of sunlight illuminating the belt, the light flux arriving 
at the ground from the belt can be calculated as follows. 

The light power scattered from a single dipole will be 


Pp=PsunX On; 


where ®sun is the power flux of sunlight and oc, is the 
optical scattering cross section of a single dipole. The 
light intensity arriving on a unit area at the earth’s 
surface from a single dipole will be 


®,= Pra/4rR?= sun XonXa/ 40k’, 


where P, is the power scattered by a single dipole, a is 
the atmospheric attenuation factor, and Ris the distance 
of the belt from the earth. The total light intensity on a 
unit area at the earth’s surface coming from a unit area 
of belt surface is given by 


Punit area of belt = Psun X on KaX DX p/A4rR?, 


where D is the depth of the belt and p is the density of 
the dipoles in the belt. 

If we assume one square degree of sky filled with 
dipoles at a density equal to that of Table III (3107 
cm’/km’), the equivalent flux is given by 


(Ve 
[57°(1/R) P 


R2 
=—unit area of belt 
Se 


PDeq deg= Punit area of belt X 


R? PsunXonXDXpXa Psun XOnXDXpXa 
Pq deg=—— X = 
An R? An (57)? 
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For the case of the test belt 


on=1.5X 10 cm?=1.5X 10-7 m’, 
D=4X10* m, 
p=20X10- per mi’. 
Thus 
Dyq deg = sun X 2.91072. 


But solar flux at the surface=®,una. The ratio of the 
two will be 2.9X10-. In other units, the equivalent 
flux at the earth’s surface is about equal to 1.6 tenth- 
magnitude stars per square degree for the test belt, 
where a tenth-magnitude star has a light flux equal to 
1.77X10- that of the sun. 

Now the effect of these light fluxes on astronomical 
observations can be estimated to be very small because 
of diffuse astronomical light sources already present. 
The average night-sky background has an intensity of 
200-300 tenth-magnitude stars per square degree. The 
zodiacal light alone ranges from 40 to 300 tenth- 
magnitude stars per square degree (Allen 1955). It may 


W. E. (MORROW TRG SNe Dar Dr Gs MacLELLAN 


be conservatively estimated that the dark parts of the 
night sky have a brightness of at least 150 tenth- 
magnitude stars per square degree. Thus scattered 
sunlight from the test belt will be extremely difficult to 
detect. In addition, at 60 days lifetime the test belt will 
occupy only a very limited portion of the sky (about 
0.15° wide as seen from the earth). Of course, it should 
be remembered that the light intensities will be much 
stronger during the initial weeks of the belt formation 
when the dipoles are closely spaced. 


RADIO ASTRONOMY EFFECTS 


Radio astronomy observations are carried out at a 
number of frequenciés ranging from wavelengths of 
several meters down to a centimeter or less. Conceivably 
dipole belts could affect these measurements in the same 
ways that optical astronomy measurements might be 
disturbed : 


(1) Attenuation of electromagnetic radiation passing 
through the belts. 

(2) Angle-of-arrival fluctuations. 

(3) Interference caused by reflection of solar or 
terrestial radio signals from the belt. 


The attenuation effects will be greatest, of course, at 
the dipole half-wave resonance since the scattering 
cross section per dipole is largest at that wavelength. 
An examination of Table III shows that the largest 
cross section occurs at resonance and is equal to about 
32 square centimeters per cubic kilometer for the case 
of the test belt 60 days after launch. Multiplying by 
the maximum thickness of the belt for this case gives 
a figure of 1.3 10~" as the fraction of 4.0-cm radiation 
attenuated in passing through the belt. 

Thus at most the 4.0-cm wavelength signals passing 
through the belt are reduced by less than 1 part in 10°. 
The total vertical attenuation through the troposphere 
at 8000 Mc is approximately 0.05 db, which is a value 


some 10‘ greater than the belt attenuation. At other . 


wavelengths the attenuation of the dipoles will, of 
course, be much smaller. 

Angle-of-arrival fluctuation effects produced by the 
belt can be expected to be extremely small and of the 
order of 10-® times the angular width of the belt as 
viewed from the ground or about 10~* deg. Fluctuations 


in signal strength produced by atmospheric variations 


can be expected to mask such small effects. - 
Interference to radio astronomy measurements can 
be caused by either natural or man-made signals being 
scattered from the belt. The strongest natural sources 
are the radiation from the sun and the earth. As in the 
case of attenuation effects, the strongest interference 
effects will be produced at the dipole half-wave reso- 
nance. The dipoles can be expected to intercept and 
reradiate electromagnetic energy radiated from sources 
such as the earth, sun, moon, and sky. The strength 
of the reradiated signals can be determined from a 
knowledge of the effective radiating temperature of 
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these bodies, together with the solid angle that they 
occupy, relative to the dipoles. In the following calcu- 
lation it is assumed that the earth and moon have an 
effective temperature of about 300°K, that the sun has 
an average effective temperature of about 10 000°K 
at 4 cm and during disturbed periods a maximum 
temperature of 10°°K. The remaining sky is assumed 
to be effectively at O°K. The angular sizes of the various 
sources are calculated for the case of a dipole 4000 km 
above the earth (see Table IV). 

Since, at resonance, the reflectance of the sky in the 
direction of the belt occupied by dipoles is only 
1.3X10-’, the increase in background temperature 
caused by the belt will be only 36X1.3X10~-7 °K or 
4.7X10-° °K on the average and occasionally may 
reach 5.2 10-6 °K. 

Such temperature changes are very small relative to 
those measurable by contemporary X-band radiometers 
which have a sensitivity of about 10-? °K. While it is 
possible that further improvements made be made in 
equipment sensitivity, variations in the atmospheric 
density probably will set a limit of about 10-* °K 
(Deem and Fanin 1955; Herbstreit and Thompson 
1955). Thus it does not appear that this source of 
interference will be noticeable, even when looking 
through the belt at its resonant frequency. At other 
frequencies the effect will be much smaller, of course. 

The man-made electromagnetic sources that can 
illuminate the dipoles are considerably stronger than 
the natural ones. These sources can be divided into 
two groups: orbital scatter ground equipment aimed 
at the belt and other microwave equipment which may 
occasionally by chance illuminate the belt. 

Typical equipment might consist of 10 kw of average 
transmitter power and 60-ft diameter antennas for the 
ground equipment. The interference to radio astronomy 
measurements will depend upon the size and aiming 
of the radio telescope and the frequency and bandwidth 
of the radiometer employed. For the purposes of this 
study, the following assumptions will be made about 


the radio telescope: 


100 ft 
8000 Mc 
10 Mc or 1000 Mc 


Antenna size 
Operating band 
Bandwidth of receiver 
Antenna gain 

Main lobe 

Far sidelobes 


65 db above isotropic 
0 db above isotropic 


The increase in the background noise temperature 
observed by a radio telescope can be calculated from 
the following equation: 


AT = PrGrWo,/(4m)2R2BK, 


where AT is the change in effective temperature, Pr 
is the average transmitted power, Gr is the transmitting 
antenna gain, \ is the wavelength, R is the distance from 
the transmitting antenna to the spot on the belt, K is 
Boltzmann’s constant, and g, is the reflectance of that 
portion of the belt illuminated by the transmitting 
antenna. 


the 


TaBLeE IV. Angular sizes of various sources for a 
dipole 4000 km above the earth. 


Angular 
aperture 
Effective (Fraction 
temperature of sphere) Product 
Earth  300°K 0.12 36°K 
Moon 300°K 4.01076 0.0024°K 
Sun 10000°K av 4.0X10-° 0.024 °K av 
10° °K max. 40 °K max. 
Sky 0°K 0.88 Oars 
Average 36.0°K 
Maximum 40.0°K 


The above value of effective temperature is valid 
for the case of a radio telescope beamwidth which is 
smaller than the angular size of the illuminated spot. 
In the case of a radio telescope with a beamwidth larger 
than the illuminated spot, the observed increase in 
effective temperature will be reduced in proportion to 
the square of the ratio of radio telescope beamwidth 
to the angular size of the spot. 

The value of o, is determined by the depth of belt 
that is illuminated and also by the relative polarization 
of the transmitting and radio telescope antennas. The 
average value of o, taking into account the above 
geometric factors is of the order of 4 10-%. 

The increases in the background noise temperature 
observed by the radio telescope have been calculated 
for various radio telescope antenna aiming conditions 
and receiver bandwidths. 

The figures in Table V are accurate to within a factor 
of 2 or 3, depending on antenna-aiming geometry, and 
indicate that in the case of the test belt the interference 
that might be experienced with radio astronomy 
measurements in the 8 kMc band because of the ground 
equipment is very small indeed. The far sidelobes 
interference is only 3X10~® °K, a level which is 3 to 4 
orders of magnitude below the present sensitivities and 
well below atmospheric limitations. The probabilities 
of main lobe interference are very low and can be 
avoided altogether since radio astronomers will have 
knowledge of the belt-orbit parameters and operating 
frequencies. 

Interference produced by ground equipment other 
than orbital scatter installations is somewhat difficult 
to estimate because of the variety of equipment that 
exists. However, some calculations can be made for 


TaBLe V. Increase in background noise temperature. Test belt 
with 10 kw, 60-ft antenna illumination. 


Effective 

Radio telescope temperature 
10 Mc bandwidth 

Main lobe O° K 

Far sidelobes 1X1078 °K 
1000 Mc bandwidth 

Main lobe 0.03°K 

Far sidelobes 1510-8 °K 
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TaBLeE VI. Radio astronomy interference from radio and radar equipment. 
Increase in sky temperature over 10 Mc bandwidth 
Main beam of 
interfering antenna: On belt Onbelt Off belt Off belt 
Radio telescope 
Belt illumination antenna: On belt Off belt Onbelt Off belt 
No. In- Average Temp. Temp. Temp.) Temp: 
Equipment stalled Frequency power Antenna OK OK OK oK 
Tropospheric 100 1000 Mc 10kw 40 db fixed 5x10-§ 1079 510-8 10-2 
communication 
Tropospheric 100 2000 Mc 10kw 40 db fixed 2K107& 10729 2GhO ee eee 
communication 
Tropospheric 100 8000 Mc 1kw 40 db fixed 2X10°* 6XK10>* 2a wo On 3 
communication 
Ground-to-air 100 300 Mc 20kw 10 db omnidirectional 3X10" 6x10 
communication Ee 
Radar search 10 400 Mc 300kw 44 db scanning 10-4 10-8 
Radar search 100 1300 Mc 10 kw _ 30 db scanning 1076 10-1 
Radar search 100 3000 Mc 10kw 30 db scanning 3X10-° Oe 
Radar tracking 100 5000 Mc 1kw 40 db tracking 31075 10-2 3x10 3107-4 
Radar tracking 100 10000 Mc 1kw 40 db tracking 10-5 10-2 1077 11078 


different typical cases. The effects produced by this 
equipment will be strongly dependent upon the 
operating frequency since the cross section of the 


dipoles varies so rapidly with frequency. Calculations: 


have been carried out on the effective increase in sky 
temperature produced by different communication and 
radar equipment which are typical of the more powerful 
presently existing or planned installations. While this 
compilation is not complete, it is believed to be fairly rep- 
resentative. In all cases a 100-ft diameter radio telescope 
antenna is assumed which has far sidelobes of isotropic 
gain. The receiver radiometer bandwidth is assumed 
to be 10 Mc. The very pessimistic assumption is made 
in each case that all the interfering equipment operates 
on the receiving frequency of the radio telescope. 
Levels of interference have been calculated for several 
different conditions of interfering and radio telescope 
antenna pointing. 

The results of these calculations, which are shown in 
Table VI, can be expected to have an accuracy of only 
about 1 part in 2 or 3 because of the various possible 
ways the two antenna beams can intersect on the belt. 
In all cases an attempt has been made to calculate 
the worst possible case. : 

Table VI shows that the interference to radio 
astronomy measurements produced by various types 
of radio and radar equipment is ‘very small indeed. 
In fact, the levels calculated are well below the best 
radiometer sensitivities achieved to date. 


INTERFERENCE TO RADIO AND RADAR SYSTEMS 


The thorough analysis of all the various types of 
interference that might possibly be produced among 
radio and radar systems by a belt of orbiting microwave 
scatterers is very complex. It is evident from the 
previous discussions that received signals significantly 
above the background noise can be produced only 
through the use of high-gain antennas and moderate 
to large amounts of transmitter power. Further, the 


strong dependence of the scattering cross section on 
frequency will tend to limit the interference effects to 
frequencies in the vicinity of the dipole resonance. 

In order to obtain an idea of the potential interference 
that an orbiting belt might produce, calculations have — 
been carried out for various typical equipments having 
high-power output and sensitive receiving systems. 
Also included is a powerful research radar. For com- 
parison the ratio of interference to receiver noise is also 
given. 

An examination of Table VII shows that the level of 
interference produced by the test belt is negligible in 
all cases except, perhaps, for extraordinarily powerful 
and sensitive research radars. 

An operational belt could employ an order of magni- 
tude greater dipole density at a higher altitude than 
that of the test belt. Even under these circumstances 
the interference levels will not increase enough to be 
observable on operational radars. However, it would 
appear that very powerful research radars could - 
observe signals reflected from the belt on those occasions 
that the belt is in the beam of these radars. 


IMPACT CHARACTERISTICS OF DIPOLES 


An important consideration is the impact hazard that 
the dipoles present to an object passing through or 
traveling in the test belt. When a straight dipole, 
moving at orbital velocity, impacts a surface normally 
(the dipole aligned along its velocity vector which is 
normal to the surface), the maximum impact effect will 
occur. In this case, the dipole will penetiass into the 
surface material a distance 


Dp=L(pa/ps)*, 


where L is the dipole length, pz the dipole density, and 
ps the surface material density. The probability that a 
dipole will be aligned closely enough with its velocity 
vector on impact to achieve this effect is extremely 
small (approximately 10-°). In fact, it is highly unlikely 
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Taste VII. Calculations for various equipment. 
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Equipment 
UHF television 


L-band radar 


S-band radar 


400 Mc radar 


C-band tracking radar 


Research radar—400 Mc 


Interference signal 
received from 


Technical characteristics test belt 


Ratio of inter- 
ference signal to 
receiver noise 


1-megawatt ERP 

10-db receiving antenna 
5 Mc bandwidth 

3-db noise figure receiver 
10-megawatt peak power 
30-db antenna gain 

100 kc bandwidth 

1 db-noise figure receiver 
10-megawatt peak power 
30-db antenna gain 

100 ke bandwidth 

1-db noise figure receiver 
10-megawatt peak power 
42-db antenna 

500 cps bandwidth 

1-db noise figure receiver 
1-megawatt peak power 
43-db antenna gain 

1 Mc bandwidth 

2-db noise figure receiver 
2.5-megawatt peak power 
60-db antenna gain 

500 cps bandwidth 

1-db noise figure receiver 


1076 w 


1072! w 


4X10! w 


10-2) w 


7X10-% w 


3X10 w 


—124 db 


56 db 


50 db 


23 db 


50 db 


8 db 


that the dipole alignment will be sufficiently close to 
its velocity vector to cause any appreciable penetration. 

Figure 4 is a plot of the likelihood of a dipole pene- 
tration of depth D, on a one square meter surface in the 
belt volume. Also shown in Fig. 4 is the likelihood of 
penetration due to micrometeorite collisions. The first 
of these micrometeorite curves is based upon micro- 
meteorite data measured by Explorer I, together with 
penetration depth information in the Proceedings of the 
Third Symposium on Hypervelocity (1959). The second 
curve was obtained from Table 27-4 in Space Tech- 


"SPACE TECHNOLOGY= 
TABLE 27-4, MODIFIED 


STRAIGHT DIPOLE 


PENETRATION DEPTH (D,)cm 
J 
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PROBABILITY OF A PENETRATION 3D, PER METER™-SEC \N THE DIPOLE BELT 


Fic, 4. Penetration depth probabilities. 


nology (1959), with the collision probabilities increased 
by a factor of 3 as suggested in the text. 

It is apparent from Fig. 4 that, at penetration depths 
less than one millimeter, the dipole collision problem is 
completely masked by the already present micro- 
meteorite collision problem. Above 1-mm penetration 
depth, an object that stays in the belt volume 
apparently faces a greater hazard from dipoles than 
from micrometeorites. In fact, however, the dipoles will 
have some curvature (see Fig. 2). For a 20 cm radius of 
curvature the maximum depth of penetration of a 
dipole is reduced to the order of one or two millimeters. 
In addition, in any practical case, a satellite, rocket, 
etc. would spend only a minute portion of its time in 
the belt volume, which effectively reduces the dipole 
hazard at least to the level of the micrometeorite 
hazard. 
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This report is divided into three parts. The first summarizes the predictions of the effects of Project West 
Ford on optical astronomy, and includes comments by this writer and by others with whom the writer has 
consulted. Part II contains several recommended investigations which might be carried out before the 
launching date which would reduce some of the uncertainties mentioned in Part I. The third section describes 
the observations which should be made of the orbiting dipoles by optical observatories the world over. 


I. PREDICTED OPTICAL EFFECTS 


OLLOWING the September 1960 URSI Assembly 

at which Morrow and MacLellan first described 
Project West Ford and the effects of the test belt on 
optical astronomy, van de Hulst and Volders in- 
dependently considered the dipole belt and published 
their findings in a report of the Inter-Union Committee 
on Frequency Allocations for Radio Astronomy and 
Space Science (October 17, 1960). The results of the 
two pairs of authors were initially in disagreement 


because, as this writer pointed out, the work of Morrow- 


and MacLellan included a factor of 10 arithmetical 
error. Changes in the West Ford experiment plans 
have made it necessary to revise the calculations of 
van de Hulst and Volders, and the new results are 
given here. 

As is now correctly stated in the accompanying 
Morrow-MacLellan report, after 60 days the visual 
surface brightness will be 1.7 tenth-magnitude stars per 
square degree. For comparison, the darkest parts of 
the sky seen through filters eliminating the brighter 
night-sky emissions have a brightness between 150 and 
200 tenth-magnitude stars per square degree. 

As van de Hulst and Volders pointed out, surface 
brightness values of this order of magnitude ‘“‘would 
not spoil now feasible photometry of faint astronomical 
sources by more than 5 or 10 percent, and that only ina 
limited area of the sky.” However, the writer must 
emphasize that these authors speak of the sky seen 
from the earth’s surface. The effect on satellite obser- 
vations will be discussed later in Sec. I. 

Since the test belt altitude remains classified, many 
important quantities cannot be calculated exactly. 
However, if one assumes the belt will be at an altitude 
of 3000 km, then at a latitude of about 40°, the apparent 
motion of the belt when in the zenith relative to the star 
background will be about 24 minutes of arc westward 
per minute of time. If one accepts the belt width at 
the end of 60 days given by Morrow and MacLellan 
(8 km), the angular width of the belt at this 3000 km 
height will be only 9’. Putting this result in a slightly 
different way, the test belt, when in the zenith at the 
end of 60 days, will appear to move its own width in 
about 23 sec. It appears that only the very fastest 
cameras using high-speed emulsions will have a chance 


of detecting the presence of the test belt 60 days after 
launching and then only by tracking the motion of the 
belt. Photoelectric devices should, however, be able to 
measure the contribution of the belt to the general sky 
background at this same time. 

The surface brightness of the test belt would increase 
near the direction of the sun, according to van de Hulst 
and Volders, because of diffraction. If the needles were 
to break up into smaller grains, as Morrow and Mac- 
Lellan have suggested might be desirable, at all places 
the belt would appear brighter, by at least an order of 
magnitude. Only the rapid rate at which solar radiation 
pressure would bring the grains down into the 
earth’s atmosphere would become, astronomically, 
advantageous. 

It has been suggested that the dipole needles be 
coated with a black substance. Even so, there might 
very well be some wavelength where the albedo ap- 
proaches unity, especially if one considers, as one must, 
the ultraviolet region of the spectrum. 

Astronomical instrumentation on board satellites 
moving in orbits under the dipole belt will, of course, 
be no freer from the influences of the needles than any 
earthbound observatory. In fact, the possibility of 
higher albedos occurring at shorter wavelengths, the 
existence of skies made far darker because of the absence 
of airglow and man-made lights, and possibly other . 
factors, will make the influence of the dipoles on satellite 
astronomy more important than on earth-bound 
astronomy. The chance of seeing for the first time 
extremely faint objects with satellite telescopes possess- 
ing resolving powers no longer limited by an atmosphere 
in a sky perfectly black except for the feeble zodiacal 
light at the ecliptic poles appears to be threatened by 
Project West Ford. The brightness of the space between 
the stars as seen from above the earth’s atmosphere in 
the direction near the ecliptic poles can only be esti- 
mated; it is probably no more than 5 tenth-magnitudes 
star per square degree. Only high-altitude satellites 
would then be unaffected by orbiting dipoles. 

As has been emphasized by van de Hulst and Volders, 
there are other serious consequences of the belt. 
Determinations of both the color and degree of polari- 
zation of astronomical radiation involve measurements 
of a very high precision. Babcock, Hiltner, and others 


114 


Wi eh OR DAS LRONOMICAL EFREECTS 


have achieved precisions in polarization measures 
amounting to better than 0.001 mag. (0.1%)... Any 
tendency of the 2-cm needles to align themselves 
parallel to each other, as may conceivably happen as 
they move along the earth’s magnetic field, could very 
well introduce serious errors in similar observations 
made in the future. That such a tendency will exist has 
been confirmed by Professor W. H. Furry of the Physics 
Department at Harvard University in conversation 
with the writer. The forces involved are difficult to 
estimate since one must know the size of the electric 
charge which each needle will undoubtedly accumulate 
through the action of the photoelectric effect and 
collisions with ions, electrons, and other charged 
particles. 

It is not possible to overemphasize the concluding 
remarks of the van de Hulst-Volders report. The 35 
kg of needles will, at best, be barely detectable with 
presently existing optical instruments. Substantially 
increased dipole densities, as the operational belts will 
very likely have, will make any research on faint 
astronomical objects most difficult except from above 
these belts. Optical astronomers must carefully evaluate 
the effects of the orbiting dipoles in the test belt as 
accurately and as quickly as possible so that an upper 
limit can be set on the size of future experiments, if any. 

Tn conclusion, it should be strongly emphasized that 
all the calculations described in this section are based 
on the values of needle density and rate of dispersion 
given by Morrow and MacLellan. It is not possible to 
check these important quantities, because the details 
of the orbit and of the dispensing mechanism remain 
classified. This is most unfortunate, since the correctness 
of the predictions made here are of vital concern to all 
optical astronomers in the world. It would be most re- 
assuring to have an independent check made of the com- 
putations of Morrowand MacLellanso that the calculated 
brightness of the dipole belt as given in their report can 
be verified. It must further be emphasized that improper 
operation of the dispenser, or the attainment of an 
orbit other than the one intended, could also lead to 
greatly different results. 


Il, RECOMMENDED PRE-LAUNCH INVESTIGATIONS 


Declassification of the details of the intended orbit, 
of the needle dispensing mechanism, and of the time 
of launching would be, without doubt, the most 
valuable action which could be taken before the launch 
date. This would enable independent calculations to be 
made of the belt density and extent, both of which bear 
vitally on the results of Part I. Even with full in- 
formation on such details, it still may not be possible to 
calculate accurately all the contributing effects. For 
example, F. L. Whipple of Harvard-Smithsonian has 
told the writer that he estimates the belt will diffuse at 
a rate faster than that calculated by Morrow and 
MacLellan. 


AS 


Albedo and diffraction measures of the various needle 
materials under consideration would be valuable, 
although the assumption that the albedo will be unity 
at some wavelength will not change the results of Part I 
by more than a factor of 2 or 3. 


TI. RECOMMENDED POST-LAUNCH INVESTIGATIONS 


Observational programs to monitor the effects of 
Project West Ford, particularly during the first few 
weeks after launch, should be established at optical 
observatories all over the world. It would be particu- 
larly desirable if the launching were to take place 
during the dark of moon, i.e., within a few days of new 
moon, preferably before. A sky unilluminated by moon- 
light would allow more precise measurements to be 
made. The writer is currently making efforts to have 
this matter taken into consideration. However, progress 
is difficult since the date of launching remains classified. 

The following paragraphs outline briefly several 
important recommendations. 

Fast, wide-angle cameras using filter-emulsion combi- 
nations which reject bright air-glow radiations such 
as \A3914, 5577, 6300, 6363, etc., should follow the 
growth of the cloud with time. All photographs should 
be calibrated either by exposing the emulsion to a 
suitable set of known intensity areas, or by photo- 
graphing under as similar conditions as possible an 
extended source of continuous radiation with an 
emulsion from the same batch as the original. Suitable 
calibration objects are the Andromeda galaxy and its 
satellites M32 and NGC 205, the Magellanic Clouds, 
and any other large galaxies. Globular clusters and 
planetary nebulae would be less desirable but acceptable 
calibration sources. 

Photoelectric photometers on land or in orbit, exposed 
to the sky ahead of the approaching belt will give 
valuable intensity data. Multicolor observations would 
be particularly useful. Care should be taken to expose 
through apertures of known size, so that the raw data, 
when combined with observations of standard stars 
and atmospheric extinction eurves, can be converted to 
convenient surface brightness data. 

Polarimeter measurements both of starlight passing 
through the belt and of the belt radiation itself are 
urgently needed. An effort certainly should be made to 
find out if there are belt orientations which produce 
more polarization than others. 

Visual observations alone would be very much 
worthwhile. Any unusual effects, such as might be 
caused by preferential scattering from aligned needles, 
should’ be reported. Naked-eye observations made 
through a sheet of polaroid rotated slowly about the 
optical axis would be particularly desirable. 

Many other types of optical observations can be 
imagined, and astronomers are encouraged to make 
them. Measurements are urgently needed, so that the 
effects of Project West Ford can be evaluated. Please 
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send all results to the writer, who will see to it that all 
measurements are carefully considered and analyzed 
and that the final evaluation is made known to all 
astronomers, as well as to the National Academy of 
Sciences and the Space Science Board. 


WILLIAM LILLER 


If the date of the West Ford launching is declassified 
as early as several days before the actual launching, 
this information will appear on a Harvard Announce- 
ment Card. Otherwise, interested astronomers should 
watch their newspapers for word of the launching. 
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The future conduct of radio astronomical research is confronted with the consequences of an orbiting 
medium surrounding the earth. Deliberate and inadvertent radio transmissions will illuminate the medium, 
presenting new and variable artificial sources of radio noise to radio telescopic observation. Basic relation- 
ships are developed permitting an estimate of the apparent brightness temperature of the medium produced 
by terrestrial radio illumination. Provided with an assumption of the power, frequency, and bandwidth of 
the transmitter, the cross section of the orbiting medium, and the properties of the observing radio telescope, 
the resultant brightness temperature may be calculated. Consequences ranging beyond the production of 
artificial radio sources have been considered. These include the attenuation of cosmic radio sources, space- 
probe radio emissions, and signal-to-noise ratios in space telecommunications systems. Based upon these 
considerations, the need for international protection of frequency bands for basic science is again evident and 


still awaits decision and action. 


d O evaluate the influence of communications 

systems of the orbiting scattering type on radio 
astronomical research, an estimate must be made of the 
brightness temperature of the medium. Consider the 
case where the medium is illuminated by a terrestrial 
transmitter, and the resulting brightness temperature 
of the illuminated region is measured by a radio 
telescope. In this discussion, antenna efficiencies will be 
ignored. Assume that the beamwidth of the transmitter 
is larger than the beamwidth of the radio telescope. 

The power per unit area at the medium is 


P.G./4nr = P,Ai/NP, (1) 


where P; is the transmitter power in watts, G; and A; 
are the gain and area of the transmitting antenna, ) is 
the wavelength of the transmissions, and 7 is the 
distance to the medium from the transmitter. 

The power per unit area at the earth resulting from 
the illuminated medium is 


(PA o/ Nr) (A mo/ dar”), (2) 


where A,, 1s the illuminated area of the medium and o 
is the cross section (area per unit area) integrated 
through the depth of the medium. 

Taking into account the collecting area of the radio 
telescope A pz, the collected power defines the equivalent 
brightness temperature: 


1 PA; Amo 
ret (CA\()sn 0 
kAve\ 7° Arr? 


where & is Boltzmann’s constant, Ave is the receiving 
bandwidth of the radio telescope (matched to the 
transmitter bandwidth Ay,), and Az is the area of the 
radio telescope antenna. In Eq. (3) it has been assumed 
that the distances to the medium from the transmitting 
antenna and the receiving radio telescope are the same. 
However, the true brightness temperature of the 
illuminated medium is fundamentally dependent upon 
the distance from the transmitter. Provided the illumi- 
nated medium fills the beam of the radio telescope, the 
brightness temperature does not depend upon the range 
of the radio telescope. 

By making further use of approximate relations © 
between Am, Qe (the solid angle of the radio telescope 
beam), and Ge (the gain of the radio telescope), we 
obtain the brightness temperature of the medium: 


Tp=(o/kAve) (P:At/4a7’). (4) 


The use of Eq. (4) requires Qe: The radio telescope 
beam is equal to or smaller than the transmitter beam. 

In radio astronomical research, it is necessary to 
distinguish, between brightness temperature and 
antenna temperature. Two differences may occur to 
reduce the brightness temperature and provide an 
antenna temperature 74. The first condition is Qr>Q; 
and the second is Avr>Av: When these conditions 
apply, then 


T sZZT 5 (Avi/ Avn)(Qi/Qe), (5) 


and in using Eq. (5) it is understood that Tz is the 
brightness temperature produced when Av,=Aver. For 


WiET Se ORRED — RAD LOM EEE CTS 


either Tz or T'4, it is assumed that the medium fills the 
transmitter beam. 


EVALUATION OF THE PROPERTIES OF THE PROJECT 
WEST FORD TEST BELT 


To evaluate the antenna temperature of the test 
belt discussed by Morrow and MacLellan (Astron. J. 66, 
107, 1961) we use the following values: P;=10‘ w; 
r=nX10° m; c=4X10-* (the average value applicable 
to the test belt quoted by Morrow and MacLellan); 
and 7,=10 m, from which 


T 4=7.2X108n= (Avg). (6) 
For a bandwidth of 10 Mc, 
Ta=72n~ °K, (7) 


where m is the altitude of the medium in units of 
thousands of kilometers. The exact altitude is not given 
by Morrow and MacLellan except that ‘‘a few thousand 
kilometers” is mentioned. If ~ is ~4 or 5, the value 
given above in (7) is in reasonable agreement with the 
3°K value given by Morrow and MacLellan in their 
example. 

To summarize the case for the test belt, the use of 
expression (6) suggests that for Ave~10 Mc, the 
illuminated belt does not pose a serious problem to radio 
astronomical research. It should be mentioned that the 
effects become more serious the more narrow the 
bandwidth. 

The dependence on bandwidth can cause serious 
effects at frequencies much lower than resonance. In 
spectral-line research, where narrow bandwidths must 
be used in many problems, the \~* dependence of « does 
not eliminate the interference hazard. This problem is 
considered below. 


POSSIBLE EFFECTS OF FUTURE ORBITING MEDIA 


To exploit the communications potential of orbiting 
media, it is reasonable to extrapolate to heavier pay- 
loads dispersed in space. In addition we may confidently 
expect larger antenna systems and more powerful 
transmitters. These advances will increase the problems 
for the future conduct of radio astronomical research. 

To estimate the influence of an operational belt, we 
assume three changes which might reasonably be 
expected in the future: 


(1) A 120-ft antenna for transmission and a 140-ft 
antenna for radio astronomical observations (both 
planned) ; 

(2) P:=20 kw (available now); 

(3) A factor of 10? increase in payload of the dipoles 
- (within the capability of planned launching vehicles). 
Using these values, expression (4) becomes 


T p=5.8X 10-2 (Ave). (8) 
~ For the test belt and a 10 Mc bandwidth we had 
Tp=72n- °K, (9) 


hy 


whereas for the operational belt assumed above we have 
Tp=58 000n~ °K. (10) 


Considering the operation of low-noise devices, an 
operational belt could be felt in sidelobes as well as in 
the main beam, creating a problem over large portions 
of the sky. It is evident that an operational belt of the 
assumed type will be of direct concern to radio 
astronomical research. 

Consider the case of the operational belt on narrow- 
band research, such as the 21-cm hydrogen line research 
efforts. 

For the assumed operational belt, we may write 


T4=5.8X 10" n-(Ave)(o21/00), (11) 


where o21/o0, the ratio of scattering cross sections, may 
be estimated from Fig. 3 of the Morrow and MacLellan 
paper. 

Since o21/o0 is approximately —40 db, or 10~, 


Dis o:5 10% * (Ave) Ke (12) 
Looking at multichannel systems where Av~6 kc, 
F497 X10 77K; (13) 


which is a substantial signal. Hydrogen-line masers are 
now in daily operation with sensitivities measured in 
fractions of one degree. Consequently much smaller and 
less powerful transmitting systems can produce measur- 
able effects. 

Consider a 10-ft L-band surveillance radar, instead 
of the 120-ft dish, and we have 


T 4 (10/120)? 9.710? nL67n-2 °K. ~— (14) 


This is based on 20 kw for a 10-ft dish. If the power is 
reduced to ~500 w, 74 remains 1.7° n~ °K, 

Again for the test belt, the spectral line problem must 
be considered. If we take the hydrogen line as an 
example, modest (60 ft) antennas and transmitters 
of ~10 kw, operating under the conditions specified 
(Avr>Ayv,), can produce 


T 4=7.2X108 n(Avr) (621/00). (15) 


For simplicity take Avrp=7.2 ke (it should be recalled 
that the multichannel 21-cm systems operate with 
banks of 5-kc filters) and the o2:/co factor of 10-4, then 


Ta=10n-22K: (16) 


which is above the threshold of detectability. 

We cannot rely on the larger beamwidth at 21 cm 
to wash out T's since much larger radio telescopes are 
under construction. It should be mentioned also that 
as \ increases, the power of available transmitters 
increases. 

Consequently, for spectral-line research at lower 
frequencies (e.g., hydrogen line, OH, deuterium, etc.), 
the \~* factor is compensated by larger antennas, 
narrow band requirements, and more powerful trans- 
mitters. The fundamental conclusion is that interference 
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to radio astronomical research can occur at wavelengths 
substantially removed from resonance. 


CONCLUSIONS 


(1) For the test belt, the estimates made above are in 
substantial agreement with the estimates given by 
Morrow and MacLellan. The effect of the test belt on 
radio astronomical research is negligible for bandwidths 
10 Mc, at resonance. 

(2) The assumed operational belt has properties 
which will create unavoidable problems for radio 
astronomical research. 

(3) For spectral-line studies, where Avg is small, 
there may be disturbances by both the test and opera- 
tional belts at resonance and at wavelengths well 
removed from resonance. The problem is not one of 
threshold detectability but one of frequency of 
disturbance. 

(4) The case for narrow-band spectral-line investi- 
gations is similar in many respects to space tele- 
communications operations where bandwidths tend to 
be narrow. As in the case of radio astronomical appli- 
cations, there will be no appreciable extinction of the 
radio signals sent to or received from space vehicles. 
However, the presence of an illuminated region in the 
telecommunications receiving beam will degrade the 
system performance and diminish the signal-to-noise 
ratio. If Tz is the receiver noise temperature, the 
system noise will be nearly 


AT~ Tr/ (Avert), 


BE) Bye iy: 


where 7 is the post detection time constant. If the 
illuminated medium is present, the system noise 
becomes 


AT~ (Tr+Ta)/(Aver)?. 


For narrow band applications and low-noise receivers, 
T can be substantially larger than Tr. If a space 
telecommunications system is accidentally operating 
through an illuminated medium, the signal-to-noise 
ratio will be degraded by 


~10 log(T4/Tr) db. 


It appears that this number can represent a significant 
loss. 

In addition to the conclusions summarized above, a 
final comment is added. The threat to the conduct of 
radio astronomical research by orbiting media should 
be viewed in perspective with other sources of man- 
made interference. New communications and navigation 
systems employing powerful transmitters in earth 
satellites will be soon forthcoming. These technological 


“advances pose another source of severe interference to 


radio research. The pursuit of basic science and the 
progress of space radio technology represent needs of 
man which must be advanced. For the impending 
interference a simple solution exists: allocation of clear 
frequency bands for basic science. This action is 
imperative and must ultimately rest on national and 
international agreements. 
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This paper contains new results on the application of the theory of periodic surfaces to the problem of 
determining the drag free orbits of satellites around an oblate earth. This theory employs a description of 
solutions which involves two distinct parts: One element describes a surface on which solutions travel and 
the second element describes the solutions on that surface. The formulae giving the approximate solutions 
were tested numerically and the results included in the paper. 


1. INTRODUCTION 


HIS paper contains new results on the application 

of the theory of periodic surfaces to the problem 

of determining the (drag free) orbits of satellites around 
an oblate earth. The basic mathematical theory, except 
for the new results of Sec. 9, is in the NAVORD report 
(Diliberto 1959a) “Periodic surfaces and satellite 


* Professors Diliberto and Kyner of the University of California 
and the University of Southern California, respectively, are 
consultants at Space Technology Laboratories, Inc. 


orbits.” The discussion given here is limited to that 
part of the theory which applies to the satellite problem. 

The proposed mathematical theory originates from a 
conjecture, the existence of a family of periodic surfaces, 
which, if true, would establish the stability of the orbits 
and the validity of the algorithm for finding approxi- 
mate solutions. This theory employs a description of 
solutions which involves two distinct parts : One element 
describes a surface on which solutions travel and the 
second element describes the solutions on that surface. 

For the convenience of the reader, this paper contains 
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an exposition of material in the report referred to includ- 
ing (i) an expansion technique for finding the periodic 
surfaces, (ii) a simple prerequisite condition which 
coordinate systems must satisfy so as to allow the 
possibility of such expansions, and (ili) the determina- 
tion of at least one coordinate system which satisfies the 
prerequisite condition. The new material includes a 
procedure for finding expansion for solutions on the 
surface (this can be described as an expansion for the 
secular terms, i.e., for the rotation number).-It was 
found that the energy integral has a simple form in 
terms of the coordinates which are used. From this it 
follows that if one component of the surface is known to 
a given order, the other component can be found to the 
same order directly from the energy integral. 

Although very little progress has been made on the 
mathematical conjecture which is the basis of this work, 
it is known that the second-order terms of the surface 
expansion are calculable (if they were not, the conjec- 
ture would be false). Once these terms and the second- 
order terms for the solutions on the surface are at hand 
(the latter are always calculable), one will have enough 
terms to make reasonable definitive comparisons with 
the various perturbation schemes which have been 
proposed. 

Despite the fact that only the first-order terms in the 
rotation number and (with the exception of certain long- 
period terms) the first-order terms in the surface 
expansions have been evaluated, it seemed advisable to 
initiate numerical studies whose eventual goal would be 
the comparison of different perturbation schemes and 
step-by-step integrations. These results are presented in 
Sec. 11. They may be regarded as experimental evidence 
supporting the mathematical conjecture. 

It is hoped that subsequent papers will include the 
derivation of the higher-order terms, a discussion of the 
qualitative properties of the orbits using these expansion 
techniques, and the longitude-time (¢—1?) expansion. 

The authors wish to acknowledge their indebtedness 
to the staff of the Mathematical Analysis Department 
of Space Technology Laboratories, Inc. and in particular 
to Dr. John Titus, who checked many formulae and 
helped in the numerical studies. 


2. THE DIFFERENTIAL EQUATIONS 


In terms of spherical coordinates (r,6,¢) the equations 
of motion of a satellite about an oblate spheroid, where 
all but the effects of gravity have been excluded, are 


dg 2 oV 
= G)-(2)-5 
or 
10V 
—r sin cosa( \- 
AER dt r 00 
1 oV 


d do 
“|r sin’6 (~) = a 
dt dt rsinO d¢ 


id 


r dt 


(2.1) 
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where V is the latitude dependent gravitational 


potential, 
GM,R R\ 71 
V= | +1( ) —cos' (2.2) 
Ry: r 3 


G is the gravitational constant, M the mass of the earth, 
R the earth’s equatorial radius, and J is a scalar which 
is a measure of the earth’s oblateness. 

Since the earth has been assumed to be symmetric 
about the polar axis, dV/d¢=0. The last of (2.1) can be 
integrated to give 


r’ sin?6(d¢/dt)= p=const. (23) 
In other words, because of the earth’s symmetry, the 
angular momentum about the polar axis is constant. If 
p is not equal to zero, then ¢ is a monotonic function of 
t and can be used as an independent variable. It is 
convenient to introduce two new dependent variables, 
U=coté, W=1/r sind. (2.4) 

Note that by this change of variables one excludes all 
orbits which pass over the poles or through the center 
of the earth. 

It is straightforward, but tedious, to calculate that 
U and W satisfy 


CU —2.UW 

eS 
d¢* (1+U?)? 

(2.5) 
ew A 5 4 
a ee al 

d¢ (1+ U?)3 (1+02)72 (1+07?)3 
where 


A=GM/p, \=JGMR?/p°. 
The preceding change of variables seems to be due to 
Blitzer, Weisfeld, and Wheelon. The fourth-order 
system (2.5) describes the motion of the satellite as a 
function of the longitude ¢. Equation (2.3) must be 
integrated if the position as a function of time is wanted. 
Equations (2.5) can be interpreted as the equations 
of a system of coupled harmonic oscillators. This leads 
one to the geometry of the four-dimensional phase 
space with variables 


w= Pie 
42=dU/dd, 


y=W 


y2= dW /do. Sen 
Although the physical meaning of the phase space 
transformations is frequently obscure, the phase space 
geometry plays an important role in most of nonlinear 
mechanics. A convenient decomposition of the phase 
space will be discussed in the next section. 
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3. PHASE SPACE VARIABLES 


The differential equations for the phase space varia- 
bles are 


dxy/do= x» 
dx, 2hey1 
F ek eee 21 
X72)? 
bs : (3.1) 
dy;/dp= y2 
dy2 


A 
SS OT ee 
dp (Atm)! 


If the earth were spherical, then \ would be zero and 
the Eqs. (3.1) simplify to 

dx1/dp=x2, dyi/do=yo, 

da2/dp=—a1, dy2/de=—yit+[A/(1+27) 8]. 


The solution functions x;(¢), yi(¢) are coupled 
through the last equation of (3.2), However one can find 
functions a;(21,%2) so that when new variables 


i=1, 2, 


5 4 
2 Se tee oe ee eee . 
: fe (i+a,2)! 


(3.2) 


Zi= yi—0; (41,42), - (3.3) 
are introduced, the unperturbed Eqs. (3.2) become 
dz;/do= 22, 


d32/do= == le 


dx1/db= x2, (3.4) 
dx2/dp=—%1, 3 


The phase space has been decomposed into the product 
of two planes. Clearly, the solution curves in the planes 
are circles. 

The derivation of the a; functions will not be given 
here. It can be verified that 


A(i+2,)? 
(1--y?+- 272") 
Axx 
ag= 
(1-+-a?-+az%) (1-22)? 


serve our purpose. 5; 
If polar coordinates are now introduced in the two 
planes, 


(3.5) 


41=71S1NO1, 21=72 Sine, 


(3.6) 


X2=711 COS61, %22=172 COSAo, 


then the unperturbed equations become 
d0;/d¢=1, dr:/do=0, i=1, 2. (3.7) 


The perturbed equations relative to this coordinate 
system are 


d0;/do= 1+ 2O;(61,62,71,72) 


3.8 
dr;/do=XR;(61,00,71,72), t= 1, MU, ( ) 


The analytic functions 0; and R; have period 27 in 6;. 
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If one writes 
Oe sind;, 


C;=cos6;, (3.9) 


Pr(u)= (Le), 
then 
Qi= 2S7P 5 (r~S1)[72S2+A P» (r1)P_1(71S1) ], 
Ri=— 2ryS1C1P5(rvS1)[r2S2+A P2(1r1)P_i(r1S1)_], 
120.= {2S1P5(rS1)[r2Se+A Po(11)P_1 (1181) J} 
X {Are Pa(ri)Pi(riS1) 
xX [— 2C1C2P_2 (ryS1) — SoS1 (1 —rP+2r?S1) ]} 
— Sof r2S22+ 2roS0A Po(11)P_1(rvS1) 
+A?P4(r1)P 2 (r1S1)} 
X{—4Ps(rSi1)+5P7(1S1)}, 
R,= {2S1P5(r1S1)[r2eSe+A Po (r1)P_1 (r1S1) }} 
X {Ar?Pa(11)Pi(riS1) 
x [— 2C1S2P_» (r1S1) +S1C2(1 —re+2r2S1) ]} 
+Cof1r22So?+ 275A P» (r1)P_1 (7151) 
+ A?P4(11)P_2(1r1S1)} 
X {—4P (7181) +5P7(71S1)}. 


(3.10) 


4, PERIODIC SURFACES 
If p: and p» are constants, the equations 


t= pisinb:, 2:=p2sinbs, 


(4.1) 


X2= pi coshi, Z2= po COSO2, 


define a two-dimensional torus in the four-dimensional | 
phase space. It is generated by the solutions to the 
unperturbed Eqs. (3.7), ie., 


A= b+, r= Ps (4.2) 


where 2; and #; are the initial values. Each point on the 
torus lies on one and only one solution curve. Further- 
more, any solution curve which touches the surface at 
one point must in fact be a curve on the surface. This 
torus is the simplest example of a periodic two surface. 

Definition: A periodic two surface of Eq. (3.8) is the 
graph of a pair of analytic functions 


1:=S;(01,92,) (4.3) 


defined for all (61,62), and for some neighborhood of 
\=0, with period 27 in @,, and such that if 


Di Si (01,02,A), (4.4) 


then the solution functions 7;(¢), 0;(¢) taking on the 
initial values 


i=1,2, 


Vi(do)=Pi, 9:(G0) = 21, (4.5) 


ri()=SiLO1(¢),42(6),A ] 
for —~<g¢<+o. 


satisfy 
(4.6) © 


The torus which has been given as an example of a 
periodic two surface depends on two parameters p; and 
po, the radii of its normal cross sections. Clearly, as pi 


| 
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and #2 tend to zero, the torus collapses onto the trivial 
solution 


(4.7) 


The two-parameter family of tori has the property of 
completeness, which will be defined at the end of this 
section. 

The solution which can be considered to be the center 
of the nested family of tori is the unique circular orbit in 
the equatorial plane corresponding to the assumed value 
of p, the component of angular momentum about the 
polar axis. In the phase space, this solution is repre- 
sented by a critical point of the system of differential 
equations, 

The central mathematical problem which is studied 


71=0,- ro=0. 


here is essentially this: Does there exist a nested com- 


plete family of periodic two surfaces when A, the oblate- 
ness parameter, is not zero? As stated in the Introduc- 
tion, this problem is still unsolved. The approximation 


scheme which is motivated by the belief that these 


surfaces do exist is the subject of this paper. 
The periodic surfaces of the type considered here are 
closely related to the concept of periodic integrals. This 


relationship is particularly important since the periodic 


integral seems to be the proper object of study for 
existential purposes. On the other hand, the periodic two 
surfaces seem to be a more useful approximation device. 
This is discussed in the next section. 

Definition: Equations (3.8) have a pair of independ- 
ent periodic integrals if there exist two analytic func- 
tions H;(61,62,71,72,A), i=1, 2, having period 27 in @,, 
which are independent, i.e., 


0H; 
Or, 


0, (4.8) 


and are integrals, i.e., for solution functions r;(¢),:(¢), 


H{01($),02(),71(),72(¢),A ]=pi=const. (4.9) 
The H; are normalized if 
H ;(01,02,71,72,0) = 1%. (4.10) 


A pair of independent periodic integrals defines a 
two-dimensional surface, the intersection of the level 
surfaces H;=p;. This surface is the periodic surface 
determined by the three parameters #1, po, and A. One 
can obtain the functions S; by solving for the 7; as func- 
tions of 6, 6: and the three parameters. If the H; are 
normalized, then so are the S;, 1.e., 


r= Si= Pr +S; (61,02,P1,)2,d)- (4.11) 


The above discussion shows that the periodic surfaces can 
be considered to be the inverses of the periodic integrals. 

The pair of periodic integrals can also be interpreted 
as defining a change of variable p;=H; such that the 
differential equations for r; are replaced by 


dp:/dp=0, i=1, 2. (4.12) 
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This leads to a convenient definition of a complete 
normalized family of periodic two surfaces. 

Definition: A normalized family of periodic two 
surfaces for Eqs. (3.8) is represented by a pair of 
analytic functions $;(6,62,p1,p2,A), with period 27 in 6, 
such that the change in variables defined by (4.11) 
results in the differential equations (4.12). The family is 
said to be complete if S;(01,62,0,0,A)=0. 

The geometry of a complete normalized family of 
periodic two surfaces is that of the example of the tori, 
i.e., the surfaces are essentially concentric tubes around 
the point defined by r;=0. As the parameters p; tend to 
zero, the tubes collapse onto the point. 

In the next section the problem of the existence of 
periodic integrals (and therefore surfaces) will be 
considered, 


5. PERIODIC INTEGRAL EXPANSIONS 


An integral to the sytem of differential equations(3.8) 
can be considered to be a solution to the partial differ- 
ential equation 


2 (0H; OH; 
SEEN O 1 Pun ai=Onun (Sen) 
a0; ar; 


a= 


Since a pair of analytic normalized integrals can be 
written 


A,=r4+> NH n'(O1,00,71,72), (5.2) 
n=1 


the single partial differential equation (5.1) may be 
replaced by the infinite set 


OH; OH; 
=—R,;, 4=1, 2, 
00, 062 
(5.3) 
OH OH 
IMs ie Ape PP 
00; 002 
where 
2 OH 0H 
PH=-> Jo +R —] (5.4) 
aed 00; 1; 


It is interesting that the Eqs. (5.1) for the periodic 
integrals are linear, while the corresponding equations 
for the periodic surface functions are nonlinear. The 
latter equations are obtained by differentiating (4.11), 
aS; 
80. 


as; 
See ere Nt TAGs) (5.5) 
1 


The nonlinearity comes from the fact that R; and O,are 
functions of S; and Sz, e.g., 


0;=0;(01, 02, prtrsi, potrS,). 


Evidently, it will be better to study the periodic 
integrals if existence is the question. On the other hand, 


(5.6) 
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the periodic surfaces have a more direct connection with 
the solutions of the original differential equations and 
can be used for approximating purposes. The linear 
differential equations which are satisfied by the func- 
tions H,,‘ are quite simple. The following elementary 
theorem about them is basic to this study. 

Theorem: A necessary and sufficient condition that 
there exist a doubly periodic solution to the differential 
equation 


(5.7) 


is that the doubly periodic continuous function g satisy 
the equation 


ox 


1 25 


Tall 


(5.8) 


Both j and g are to have period 27 in 4; and 6». 
From this theorem it follows that the first equations 
of (5.3) will have periodic solutions if and only if 


MR;=0, i=1,2. 


Let us suppose that this equation is satisfied. Then the 
solution to the differential equation (5.3) can be 
written as 


Ay'=QR;+hi' (2.— 61), (5.10) 
where the operator Q is defined by 
1 Qa 
Qg=— iL tg(0i+t, 02+#)dt. (5.11) 
Qr 0 


The hx’ are periodic functions which are yet to be 
determined. 


Now 
0H,' dH! 
ty = POR: Phy. (5.12) 
06; O62 
The equations , 
MP(QR;+h;')=0 (5.13) 


must be satisfied if H2‘ are to be periodic. Therefore, 
substituting in (5.13) the expression for P, given in 
(5.4), one finds for /;' the differential equation 


diy’  MPOR; 


—_=—————,, where y=02—61. 
dy M(®.—®,) 


(5.14) 


Fortunately, for the satellite problem, the denominator 
is nonzero. If the iy‘ are to be periodic, then 


2x VPOR: 
0=)'(2n)—hy#(0)= Tf Pee ah pet) 
0 M (O.— @1) 
This is another condition to be satisfied by differential 
equations. By continuing in this way, one obtains an 
infinite set of conditions. 
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The functions hm’ are determined to within an | 
additive constant by this procedure. It is convenient to | 
normalize them by requiring that 1 


hm’ (ar/2)=0. (5.16) I 


6. TRANSFORMATION OF THE S2 EQUATION 


A necessary condition that the satellite equations jj 
have a complete family of periodic two surfaces is that 


Wi; =O7 =, (6.1) 
It can be shown that 
MR,=0 
MR2= [reK (r1) sin2y ]/2r 
MOQ,=11(1r1)/20 
M0O.= [Le (n)+K (rn) cos2y ]/2r, 


(6.2) 


where 
y= O2.— 61 


Ly (11) = 20 A (1+72)-? 


TA 
L2(n)= ee (6.3) 


TA 
K(n)= ins 4 (4—r2)r?. 


It is therefore necessary to seek another coordinate 
system in which an expansion is possible. 
Let us make a transformation of the 7; equations 


(6.4) 


i Si 


To 8. (01, 92,51,S2) 
so that 
ds;/db=XT ;(01,62,51,52) (6.5) 
with 


MT;=0. (6.6) 


It follows from differentiating (6.4) that 
ARi=AT1 

Og og 
AR2=—(1+AO1) +—(1+A®2) 

00; 002 (6.7) 
0g 0g 

+—ATi+—aT2. 
OS1 OS, 
The zero-order terms in will vanish if 
&=8(82— 81, 51, 2). 

Condition (6.6) will be satisfied if 


MR=M| teed, (6.8) 


| 
| 
j 
| 
| 
| 
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Since Mg=g, and M(0.— 


8 (O2— 
is taken as the solution of the differential equation 


0,)+0, the function 


1, Si, So) 


dg MR, —K(ry) sin2y as 
SSS ee g, (6.9 
dy M(@,—®0;) Li(r1)—L2(r1)—K(n1) cos2y 
such that 
g(0,51,52) = So. (6.10) 
Hence 
ro= g(02—61, $1, S2)=S2[1+s7 sin}. (6.11) 
The conditions MT;=0 are now satisfied. 
The differential equations are now 
d6;/do= 1+rA0,[ 61, Oo, Si, £(A2— 61, Si, 52) | (6.12) 


ds;/do= AT; (41,02,51,52), 


- where 


T1= Ri[61,02,81,8 (82— 151,52) | 
i [1+s? sin? (@2— 61) | 
Or» 


a epee 5 | not (=) 


It will be shown in the next section that sz is pro- 
portional to the eccentricity of the instantaneous 
ellipse associated with the orbit. 


7. THE ENERGY INTEGRAL 


The formal expansions of the periodic integrals which 
were developed in Sec. 5 seem to be independent of the 
periodic integral which is known to exist, namely the 
energy integral. In this section, it is shown that Ho, 
the integral associated with the variable so, is essentially 
the energy integral. This fact gives us a convenient 
formula for computing Ss». 

The equation r2= g can also be written 


$2 =re[ 111? sin? (02— 61) J (7.1) 
or 
= (2-++2.?)+ (a12— 221)”, (7.2) 
where, as in Sec. 3, 
w1=n sing, 
X2=11 COS61, 
r 13 
21> 172 sinfo= yi1—a1(%1,%2), ( ) 
Z2= 12 COSA2= Yyo—2 (41,42). 
Therefore 
$= yt yo? + (x1y2— a1)? +ar+ae?+p* 
—2{y1(a1+228)+y2(e2—418)}, (7.4) 
where - 
B=ay%2—O2%1. (7.5) 


If the square of the velocity is written in spherical 
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coordinates, then one has 


[HZ Q)] eo 


or 
V?= pL (a1ye— aay)’ + yr+y2" |. (7.7) 
Hence 
V2 
a ea ee 
—2{y1(ait«8)+y2(a2—218)}. (7.8) 
Using (3.5), the defining equations for a;, one has 
a2—418=0 
A 
O00 ee 
(1-+-«;7)? (7.9) 
Aay 
artar+e= 5 
GEES) 
Therefore 
V2 = 2A A? 
52=—— ae : (7.10) 


Po aGsea)e esi) 


On the other hand, the energy integral can be written 


ZEA = DAS INL 
= —— ~cos] (7.11) 
Pp Pp Ter 
Noting that 
yiA/(1+a?)!=A/r, (7.12) 
one obtains the desired result, 
A? 2h OAL 
2= Sse [=—cow]. (7.13) 
(1+s?) 2£ 7'L3 
The angular momentum per unit mass is 
dd\? db\? 
-r| (=) +-sino(—) (7.14) 
dt dt 
or 
H?=7[1+-s7]. (7.15) 
Hence 
51= tant, (7.16) 


where 7 is the inclination of the orbit plane. 
The instantaneous eccentricity satisfies the equation 


e’=1+[2E,.H?/(GM)* ], (7.17) 
where 
V2? GM ppl 
= E+ | —— cos} (7.18) 
2 r r? 13 


E is the total energy per unit mass. 
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It follows that 
5o2= e®A?/(1+5,2). (7.19) 


It has been shown that the variables s2 and s; are related 
to the instantaneous eccentricity and the inclination of 
the orbit plane. 


8. ANGULAR MOMENTUM 


The variable s; is related to H, the angular momen- 
tum per unit mass, through the equation 


B= p(1+s7). (8.1) 


A method for approximating the variable s; is therefore 
of independent interest. However, its role in this paper 
is to give us the first-order terms in the expansion of the 
periodic two surface. These will be used in the approxi- 
mation theory of the following sections. 

It is convenient to use the notation of Sec. 2. Recall 


that 
aU? 
= u+(—) : (8.2) 
do 


If the first equation of (2.5) is multiplied by dU/d¢, and 


the second by dW/d¢, then they can be written as 
1 dp sdUy? 

Pcie 

2 dol \ do 


1 
=-JRA w—|————| 
3 dol (14+-U2)! 


2 d AW 20 RAAY dW 
=—J R? | es (8.3) 
3 dplL(1+U*)?1 3 (14+U*)? do 
1 dpsdW? A dw 
AG) bate 
2 dpi \ do (1+U?)3 if 
W (1—4U?) 
pee ee (1 4. U?) 7/2 
Hence 5 
2 dW 2AW 
iG) eee Go ee 
de (+09), 
4 (1—4U2) dW 
=—-JP?R‘AW?—____ —_._ (8.4) 
3 (1-+-U?)72 dd 


After integrating (8.4) 


a) Ge 


2AW 
JR] ——— 
3° La+o 


Q+03)! 


HG) 


+We— (>) —1:]+00%. (8.5) 
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Zero subscript or superscript, e.g., (dU/dd)o or s;°, |} 
indicates initial value. Hence 


(8.6) 


s°= (s1°)?+2ALG (61,2,51,52) — Go], 
where ; 
ip 2W (dW/db)-+W? ‘ 
e=- | | (8.7) 
3L(1+U?)3 A 


We have from (8.6) 


r | 
s— si +—[G(61,62,51,52) re Go |+0(Q2), (8.8) i 


Si 
and the assumed expansion 
S1= Pi tAS1'(01,62,p1,P2) +0 (d”). 


These two expressions are related by the fact that for 
all 1, 62, 


(8.9) 


G(61,02,p1,p2) — prS1' (01,02, p1,P2) = const. 
This can be shown by substituting (8.9) into (8.8). Then | 
51= s+ (A/51°)[G(01,02,1,h2)-Go]+O(’). (8.11) 


Hence 


(8.10) 


O= pi-—si°+ (A/s1°)Go 
+ A/s1)[sYSt—-G]+0(2). (8.12) 
Since both Si! and G are periodic in 6; and 6, 
soSp= ys Qnm expLi(76i+ m2) |, (8.13) 
G=D gam expli(n6itn62) |. 
Therefore 
Gnm=£nm if n?-+m?X0. (8.14) 


Formulae (8.6) and (7.13) can be used to compute s1 
and s2 if the solutions to the unperturbed equations are 
used in the expressions which are multiplied by X. 


9. APPROXIMATION OF THE THETA SOLUTIONS 
The periodic surfaces will now be approximated to 
first order by 


Si= SPF Gi (61,92,51°,52°) —G}, (9.1) 


where the functions G; are obtained from equations 
(7.13) and (8.6). To first order, the differential equations 
on the surface are 


d0,/dp=1+2O 11, 2, 51°, g(O2— 1, 81°, 52°) J. 


Approximate solutions will be found by the formal pro- 
cedure given in a later report (Diliberto 1959b). This 
differs from the earlier procedure (Diliberto 1959a), and 
as the numerical tests show, the new approximations are 
superior. 

It will now be shown that Eqs. (9.2) can be trans- 


(9.2) 
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formed into a trivial system 


dqi o 
=1- > v\".™, 
db n=1 


(9.3) 


by a formal change of variable. The desired form of the 
transformation is 


gi= 9; +6," (02— 41) 


“Us 3 A"{Di"(62— 61) +B," (01,02)}, (9.4) 


n=1 


where the 6;” are analytic with period 27 in 6.—61, and 
the B;" are analytic with period 27 in 6; and 6. Note 
that for A=0, the differential equations are identical, 
but (9.4) does not reduce to the identity transformation. 
Once the q;(¢) are known from (9.3), the 6;(¢) are 
found by inverting Eq. (9.4). If (9.4) is differentiated, 
then one has 


n 


wo % b; 
14+>5 A%,M= 1+d0;+ALO.— 0, ] De \” 


n=1 
w OB; OB," 
(2) 
n=1 06; 005 


+r(6, +0, ) |. (9.5) 
06, 0602 


Therefore the B;" must satisfy the partial differential 
equations 


0B; OB? db? 
=c;)—[@.—@1 }—+ 9; 
00; 90. dy 
OB; OB," db;" 
=c;5%—4[O.— O01] (9.6) 
00, 062 dy 
0B; 0B; 
+0; +0, ; n= De 
06, 0602 


To ensure that the B;” are periodic, the 5,” are selected 
so that 


db 
co =—M[02— 01]+ MO,, 
dy 


(9.7) 


@Be AB 
+u| 0; +0» | n>2. 
a6 36. 


_ And to ensure that the 0," are periodic, the constants 


bel OND dl Ges Wei bane BS 


c;™ are selected so that 


27 Fo, — MO, |d 
b(n) —o00)= f ean 
0 M[O.—®,] 


ye (27) —6;7 (0) 


i 
0 


OB," 
cs Me| C+ 6, 
00; 


M[O.— 0; ] 


—] 
062 


Since above conditions are applied to the derivatives of 
the 6;", some type of normalization is required. It is 
convenient to set 0;"(7/2)=0. 

The B,;” are given by 


db,° 
Bi=0| c+ —[O.— O01} oa a, 
dy 
Ne. 
Be=0)a—[0:~ 61 (9.9) 
dy 
0B; 0B," 
-[e, +O. |}: n> 2. 
00, 082 
Using the formulae of Sec. 6, one has 
M®,=11(s1°)/2r=const. (9.10) 
Therefore 
1) = 11 (s1°)/2z, (9.11) 
6°=0. 
Let 
Ad= 62.—61, 
AB= B,!— By, 
ee ean (9.12) 
Aq= QQ. 
It follows from (9.7) and (9.9) that 
B= __0(uT@,-@,)}-[@,-.} 
AB= 2 hie eS 2—-1]5, 
M[O2— 0, ] (0.13) 


Bi= QO{MO,— Qy}. 


For the numerical tests, the expansion (9.4) was 
replaced by 


gi=O1+ABy, 


14 
Ag=A6+b,(A8)+AAB. a 


[The functions 5;1(A@) should be included here, but they 
have not yet been evaluated. ] 
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Té can be verified that 
A8+6.9(A8) 


ry 
=—ttan— HC 1+(s{F 72 ee se") ; (9.15) 
2 
The easily obtained equation 
A [4— (s3° | 
fe ee (9.16) 


2 [1+(s°?F 


has been used in the evaluation of (9.15). The functions 
8(é) are obtained by iteration by means of 


cs 
A§=——tan™ 
7 


{r+ si? } 


= (9.17) 
xtan( = as) +)AB(6;,A6) } 


#,= g:—h Bi (41,A8). 
The initial approximation is the solution to the un- 
perturbed problem, Le., #;=02+ (6—¢"). 

The details of the computation of the functions By 
and AB will not be given. However, it should be noted 
that the functions 0;— M0; can be written as thesums 
of products of the form &(Aé@)/(@:), where Mf=0. 
Therefore the formula 


& 1 2a = 
O7(4)= {. f(s\ds—— if f f(s)ds di (9.18) 
6 2x a i) 


can be used. The formulae for B;' and AB are in Sec.10. 

The approximations to the periodic surfaces (which 
are assumed to exist) and the approximations to the 
motion on the suriaces have been derived. The two sets 
of formulae give the approximate position of the satellite 
as a function of ¢, the longitude. The following sections 
contain these formulae and the results of the numerical 
tests. 


10. FORMULAE 


An IBM 709 Fortran Monitor program based on the 
formulae of this section was written. The following nota- 
tion of Sec. 3 is used: 

S;=sin6;, = t 2 
oo cos4;, — A* 2 
Pi(%)= (+2 if 
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Input of initial values and constants includes go, Xo, — 
Vo, 20, Zo, Yo, Zo, ‘fe G, M, R. Additional initial values and : | 
constants are obtained using these formulae: 


r= xe ye +20 


p=xoYo—Toyo 
A=GU/P 
A=ART 


Uo=20/ (xe? +02)? 
We=1/(ae-+ye)} 

(dU /dd)o= i= Uotxetot-yoyo) +20(x0?-+yo?)* /p 
Xoetot YoYo 


— 

=) ate * p(xP+y0?)? 
s’=[Ue+ (dU/d¢)? } 
ki=Wo—AP2(s°)P_1(Uo) 


n= (=), -40,(— —) Pas °) Ps(Uo) 
oe—tan| rs (=) } 0<6,9< 2 


6s =tan— 1(Ri/ke), 0<6°<29r 
go= (AY +k)? 
hives 3: 


“Aa 


Ly(s°)= eae we 


b= 


(BZ), and AB, are obtained by evaluating the expres- 
sions below for By and AB with 64,°, 42°, s1°, and go as 
arguments. 

For desired position of the satellite the longitude, ¢, 
mod 2x, and the number of revolutions, g, preceding the 
current revolution are used in the following manner to 
obtain rectangular coordinates and velocity components. 

Values for a spherical earth are obtained by means of 
this group of formulae: 


1=$+O9— $0 
62=$+6/— $0 
AG= 62— 4; 
U=siSi 
dU /db=s!Ci 
W=goS2+A P2(s1°)Pa(U) 


dU 
Rerun Ue ee 
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Values for an oblate earth are then obtained by AB= P2(s; sinA@) P_1(s1) 
employing the next group of formulae, in succession, er tay part ae 
using the most recently obtained value of each variable XLO{MO2— O2}—O{M8:— 833] 
beginning with those for a spherical earth. A second 6:= 0°+[1+A P;(s1°) ](6—¢0) 
iteration of the computations for this group,employing 
values obtained by the first iteration, generally increases +AL(Br)o— Br ]+ALi(s1)q 


the accuracy of the results. Changes produced by a 1 zi 
. . : A = . a 0 

third iteration are very small and apparently smaller Coca tee 

than the error due to the omission of higher-order terms Xtan{tan[P;(s,°) tan(A@—4z) ] 

in the expansion. 


—d(AB— AB o)—4AdP3(s1’) 
BY=Q{M0,— 03} 


X[4— (sx)? ](@—$0+-279)}) 


O{M0:— 3} 62= 6: +8 
SiCi[P_-1(s1)—-1] E 
casera) s2= coren{- -WP(U firs (— ) ] 
+2eC1P1i(s1S1)Pa(s1) cosA0 4 3 aw 
—29CP3(s1S1)Ps(s1) cosA0 As (Wa)+—| Wa Cy 
—2¢S3P3(s1S1) smAO+AS1C1P2(s1S1) Pa(si) U=siS1 
Q{M@.— 0} dU /db=s:Cy 
= —C P;(s1S1) cosA@LA?P10(s1) (1—s2— 2514) /¢ r= P_,(U)/W 
+gPo(s1)(1—4s1?) cos*Ae] cos’@= U2P2(U) 


—ACP»(s183)Pe(s1) (1-422) sin2A@ 
— (7/3g)A*sPC3P3(s1S1)Ps(s1) cosA@ 
—4gCP3(s1S1)P6(s1) cosAd g=[se2Po(s sind) }! 
X[-149s2-+2(2—9s2—6sx4) sin2A8] W = pSv+AP-a(s15,)Palss) 
—$AsPCYPa(spS1)Po(s1) sin2A6 
+gsC1°P5(s1S1) cosA@ 


E JR 
$2= A?P2(s1)+2A J+ (—- cos) 
GM # \3 


dW /db= gC2+ASiCis?P3(s1S1) Po(s1). 


After completing the desired number of iterations 


X[Poe(s1) cos*Ad—3Pa(s1) sin?Aé ] using the above group of formulae, the following 
: : ; formulae are employed to obtain position and velocity 
2 2 
eae emeele sitAGt A*Pi(sy)/e | in rectangular coordinates and elements of the ellipse: 
+AsPS?Po(s1S1)P4(s1) sin2A6@ 3 
x=cosd/ VW 
—45P3(spS1) sinAé 
; =sing/W 
X[7A2%s2P(s1)/g—3¢+¢(4+6s7) sin?Aé ] a 
—tAs?PSy 4P4(s1S1) Po(s1) sin2Aé a / 
— gsPS1°P5(s1S1) sinAé[3— (4+s,°) sin?’Aé ] sr 
$AC1S1P2(51S1) Po(s1) t= -a(xt/n w) 
X[—2+9s2+ (4—8s2—s14) sin?A@] 
—3As2S:C8P (5151) Po(s)[1— (2+52) sin’Ae] s 2 edd o ) 
+4A P7(si)[2—3s2+ (452—s14) sin?Aé@ ] 
SiCi[ P_1(si1)—-1 dU dw 
xeer[ i, Wz ah W ) 
: C2+S2PP_1(s1) 


us [5+4P_1(s1)—6P_2(s1)— P_s(s1) ]Pe(s1) sin2A@ == P_1(U)/W=length of position vector 


20 P_-1(s1) +1] i=tan~s,=inclination of the orbit plane 
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Taste I. Summary of numerical results. 


Comparison method 


Theory Integration E 
Encke Adams-Moulton 1077 
Encke Adams-Moulton 10-7 
Encke Adams-Moulton 1077 
Encke Adams-Moulton 10-7 
Encke Adams-Moulton 8-107 
Encke Gauss-Jackson 31050 
Cowell Milne, double 10-4 

precision 
Encke Adams-Moulton 1077 
Encke Adams-Moulton 10-7 


Max 
e a n Ar 
0.2 ils) 6 500 
3 0.19 1.36 50 15 900 
0.2 TS 3 800 
0.2 125 1 150 
20 4500 
64 14900 - 
0.2 1:5 64 17 300 
0.2 1.5 64 14 400 
0.2 ee 1 150 
20 4500 © 
: 64 14 400 
OFZ 5 3 700 
0.2 15: 6 2500 


wh 


DPE Sy) ee seal 
OE PAC at 
A GU REINS 

=eccentricity 


Ee eR fil a1 
a=— | —+ (_-cov) | =semimajor axis 
GM 3 


height at perigee=a(1—e)—R 


a 


height at apogee=a(1+e)—R. 


The remainder of the program consists of comparisons 
of the above values with those of the reference trajectory 
from which longitudes were obtained for desired times. 

Formulae of this section were derived independently 
by two people. Results were then compared and dis- 
crepancies were eliminated. 


11. NUMERICAL RESULTS 


Major numerical results for the program described in 
Sec. 10 are summarized in Table I. The maximum Aris 


given in feet and is the square root of the sum of the 
squares of the differences between the Diliberto coordi- 
nates and those of the step-by-step integration run from 
which input longitudes were obtained. Other headings 
are H, error criterion for the step-by-step integration ; 


- 7, inclination in degrees; e, eccentricity; a, semimajor 


axis in earth radii; and m, number of revolutions. The 
period of the 34.3° inclination orbit is 134 minutes. All 
other periods are 155 minutes. 

Part of the attempt of the Mathematical Analysis 
Department of Space Technology Laboratories to 
obtain accurate trajectories is reflected here by the 
inclusion of comparisons with trajectories obtained by 
more than one step-by-step method for the 45° case. 
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On the basis of a recent theory of this author, satellite orbits with very small eccentricities have been 
investigated. It was found that such orbits can be represented by the formula 


r= dx (1—e cosv+4e? cos2v-+1412s? cos2w) +0 (vr), 


where w= (1+6)v-+w. Graphical illustrations of this result are given. 


OR several practical purposes it is advantageous 
to have satellite orbits with very small eccentrici- 
ties. In the cases of satellites 1959 gamma (Discoverer 
II) and 1960 82 (Tiros I) the eccentricities were 0.008 
and 0.003, respectively, and among the pre-launch 
nominal elements we even found the value e=0. Un- 
fortunately, most of the satellite theories developed 
during recent years exhibit some difficulties if the 
eccentricity is a very small quantity, because it appears 
as a divisor in the formulae expressing the perturbations 
in the osculating argument of perigee w, and mean 
anomaly M. See, in particular, the papers by Brouwer 
(1959), Garfinkel (1959), and Kozai (1959). The diffi- 
culties, however, can be easily removed from these 
theories, as pointed out by the authors themselves. In 
the application of Hansen’s method to artificial satellites 
(Musen 1959), on the other hand, it is difficult to see the 
meaning of the otherwise so very useful approximation 
by a rotating ellipse, if the eccentricity becomes ex- 
tremely small. As a matter of fact, even the geometrical 
meaning of a very small eccentricity is not quite clear, 
because the perturbations play an important role in its 
definition as a mean orbital element. The determination 
of satellite orbits also presents some difficulties in the 
case of very small eccentricities. 

One would expect more mathematical simplicity, 
rather than more complexity, for nearly circular orbits, 
so that the difficulties just mentioned must be of an 
artificial nature. Indeed, the problem of such orbits is 
not a new one in celestial mechanics, and there are 
several possible ways to avoid the small divisor e. For 
instance, instead of using the osculating elements e and 
w, introduce their combinations, e sinw and e cosw, and 
compute the perturbations in M+, or v+-w, where v is 
the true anomaly. If one prefers to work with canonical 
equations, Poincaré’s canonical elements 


(ua)* M+wo 
{2(ua)*[1— (1—e)! ]}} sinw {2(ua)*[1— (1—e?)? J} cos 


can serve the same purpose. Or determine the perturba- 
tions in the elements a, e, w, and M separately, and com- 
bine them into expressions of perturbations in the radius 
-vector r, and in the argument of latitude »+w; the small 
divisor e will cancel out, as indicated in the papers of 
Brouwer, Garfinkel, and Kozai. Because this is so, one 


would naturally think of developing a perturbation 
theory using coordinates directly instead of orbital 
elements, as was done, for instance, by King-Hele 
(1958), Strubel (1960), and Izsak (1960). Strubel paid 
special attention to the problem of small eccentricities 
(and also to that of the critical inclination), and the 
results to be presented herve are partially contained in 
his work. His views on the merits of classical methods in 
celestial mechanics, however, are somewhat too pessi- 
mistic, in my opinion. 

The point is the following. Osculating elements are 
introduced in celestial mechanics because one would 
expect them to have much slower variations than the 
coordinates, and because the differential equations 
describing their variations are usually easier to deal 
with. In the case of very small eccentricities, however, 
this is no longer true. But, in any theory using orbital 
elements we are obliged to introduce certain mean 
elements, and the question is: How well adapted are 
these mean orbital elements to describe the disturbed 
motion? 

The paper by Izsak (1960) contains a complete 
second-order solution of Vinti’s (1959) dynamical 
problem; for our present purpose, however, we need 
consider only the first-order part of the relevant expres- 
sions. Let us introduce the oblate-spheroidal coordinates 
p, g, and a, with a linear eccentricity c, so that the 
rectangular coordinates become 


x= (p?-+c?)?(1—«*)} cosa, 
y= (e+e)! (1—o*)} sina, 
2= po, 


a being the right ascension of the satellite. The constant 
c is defined by the relation 


c=J5an, 


where J» is the coefficient of the second-degree Legendre 
polynomial in the earth’s gravitational potential, 


Kb ao ar\” z 
va En) (2) ] 
r n=2 r r 
and dz is the equatorial radius of the earth. Also, for 


the radius vector of the satellite, we have the expression 


r=[p?+e(1—0?) }}. (1) 
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The semimajor axis a, the eccentricity e, and the sine 
of the inclination s=sinJ were defined as follows. In 
the exact solution of Vinti’s dynamical problem the 
coordinates p and o turn out to be even and odd 
periodic functions of the single variables » and w, 
respectively; therefore their minimum and maximum 
values, 

Pi, p2 of p 
and 

—o1, 010i o 


will remain the same during all revolutions. Let us 
define 


P2— Pi 
a=(pitp2), e= , and s=o1. 


pit pe 


The argument v is analogous to the true anomaly in a 
Keplerian motion, and w=(1+e)v-+w corresponds to 
the argument of latitude; here the quantities « and w are 
constants, the first denoting the motion of perigee. We 
shall also use the abbreviation 


v=c/4a, 


and an eccentric anomaly £. Then the first-order 
solution of Vinti’s dynamical problem can be written in 
the simple form 


p=a(1—e cosk) (2) 


DHE, 
o=Ss su(—wir') (3) 


INS, 
a—Q.= tan] ds")! n(—wl2) | 
Tv 


3°(1—s*)# — 2v2(1—52)3 


——$——— € sinv 
2(1—e?)2 (1-—e)? 
v(1—s?)3 
————- sin2v+O(v4) * (4) 
4(1—e)? 
v?(1—s?) 
M.=A(l—1.)=B-|1- | sin 
1-=¢ 
p52 ‘ 
——— sin2w+O(r‘), (5) 
4(1—e)3 


where (2, and ¢, are constants of integration, MV, and # 
denote the mean anomaly and the anomalistic mean 
motion of the satellite, and sn[(2Z/m)w|??] and 
tn (2L/m)w|2?] stand for Jacobian elliptic functions 
with parameters /?, The quarter real period of the func- 
tion sn(z|/) is 


mp sly? 1-3)? 
Pb Ne 2-4 
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so that the period of 


2L Bs Ie 
su(—w | r) = (1+) sinw-+— sin3w+O (v4) 
T 16 16 


becomes 27. The connection between the true and |i: 


eccentric anomalies is given by the formulae 


2K cosH—e,, 
cn (— | #) = ————_., 
T 1—e, cos# 
6) 
2K (1—e,2)} sinE 
sn{ —»| #) =———_—__, 
Tv 1—e, cosE 


where ¢x is another eccentricity, namely 


2 


ene 


Vv 
1-é 


1-235+06 | (7) 


Finally, the developments of the quantities ¢, #, k?, and 
P are, up to O(y4): 


v 3y4 
e=——_(4— 5s?) +—_____ 
A4(1—e)? 64(1—e?)4 
X[(96—432s?+345s*) — (48-+96s?—170s4)e ] 


3A(1—s*) 3yt(1—s?) 


2(1-&) 8(1—#)3 
XL(-+11s)— se] 
ye? vie? 
“Woe G2ab 
prs? 4y4s? 
a= ae (i— Ne 
je (et 


the motion of the node was found to be 


32 (1—s?)? 3y4(1—s?)} 
q= —-————_+——_ (18 — 13s”) +24 se]. 
2(1—e)?  16(1—e*)4 

In all these formulae given above the eccentricity 
never appears as a divisor, even in the mean anomaly 
(because Mx is not the osculating mean anomaly, but 
just a linear function of time). Hence they are valid for 
arbitrarily small values of e, the smallest possible value 
being e=0. Because the eccentricity was defined in a 
geometrical way, it is easy to see the meaning of the 
condition e=0, that is, p(t)=a. The radius vector 7, on 
the other hand, whose behavior determines the con- 
secutive perigees and apogees of the satellite, is by no 
means constant; rather, we have according to Eqs. (1), 
(2), and (3) the expression 


r=a| +s] 1-¢ (cir) |}. (8) 
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This formula already shows the fact that the orbit has 
two apogees with a distance of 


Tmax= a(1-+v*)* 


at the ascending (w=0) and descending nodes (w=7), 
and two perigees with a distance of 


Ymin= [1+ v?(1—s?) } 


at the highest (w=7/2) and lowest points of the orbit 
above the equator (w=37/2). (Only for equatorial 
orbits, that is, for s=0, does e=0 mean a circular solu- 
tion.) Such an orbit is very nearly an ellipse with its 
center (not a focus) at the center of the earth, having a 
semimajor axis of ¢)>=a(1+ 4’) in the direction of the 
line of nodes, and an eccentricity of es>=vs. Note that 
this ellipse does not rotate in the orbital plane. 

It is more difficult to visualize the nature of orbits 
with very small, but not vanishing eccentricities. Let us 
consider again Eqs. (1), (2), and (3), from which, 
neglecting quantities of O(v*), we obtain first the 
expression 


py? 1—s? sin?’w 
rma( 1—< cosH-+— ——). 
2 1—ecosE 


Furthermore, if we assume that e=O(v), then Eqs. (6) 
and (7) give 


e cosH=4e?+e cosy—4e? cos2v+O(r), 
so that we can write 
r=dx(1—e cosv+4e? cos2v-+4r*s? cos2w)+O(r*), (9) 


where dx=a[L1—4e?+437?(1—45?) ] and w= (1+e)v+w. 
There is no simple geometrical meaning that can be 
attached to this equation, as shown by numerical 
examples such as those illustrated in Fig. 1. The several 
curves in this diagram represent the deviations of the 
radius vector from its mean value ax, for different values 
of e and w, It should be noted that w is not the argument 
of perigee in the usual sense, but only a constant of 


‘integration in the theory; the linear function o=e«v+w 


can considered to be the mean argument of perigee. The 
reason for varying w is the following. As the true 
anomaly v increases by the amount 27/e, the phase 
difference w—v=@ takes all values between 0 and 27, 
and the shape of the orbit changes from revolution to 
revolution, except for e=0, or s=0. These changes can 
be simulated by letting v increase from 0 to 27 only, but 
varying the constant w. 

In the construction of Fig. 1 the following constants 
were used : 


a=7100 km, 
T=48°39, 
- Jy=1,0822 10-5, 


d= 7102 km, 
s=sin[ =0,7477, 
v=0,000873, 4v?s?=0,000122. 


131 


° 45° 90° (35° 180° 225° 270° 315° 360° 


Fic. 1. The variations of the radius vector as functions of the true 
anomaly, for different values of e and w. 
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In the expressions of the argument of perigee w and the mean anomaly M, derived in a previous paper on 
the motion of a close earth satellite (Kozai 1959), there is a divisor e, the eccentricity. If the eccentricity is 
very small, the amplitudes in these expressions become so large that they cannot be regarded as perturba- 
tions. The present paper proves that when the three elements e cosw, —e sinw, and M/-+w are used instead 
of e, w, and M, the expressions of the radius and the argument of latitude are the same for beth the typical 


and the small eccentricities. 


VOLUME 66, 


NUMBER 3 


——— 


INTRODUCTION 


F the eccentricity e of an artificial satellite is a small 
quantity of the first order or less, it seems that the 
results derived in a previous paper (Kozai 1959) do not 
hold because of a small divisor e appearing in the expres- 
sions of the mean anomaly M, and the argument of 
perigee, w. In this case the amplitudes of the inequalities 
having the divisor e must not be regarded as small 
quantities, but are sometimes so large that the usual 
method of successive approximation seems to lead false- 
solutions. 

In fact, as the disturbing function was expanded into 
a power series of inequalities of the orbital elements in 
the previous paper, this series may not be convergent 
for the small eccentricity even in a rough sense of 
convergency. 

To derive the mean value of the eccentricity from 
observations the mean value of the argument of peri- 
gee must be known, and to derive the mean value of 
the argument of perigee that of the eccentricity is 
necessary. For the small eccentricity it is difficult to 
know the two mean values simultaneously by using the 
previous expressions. However, after subtracting the 
short periodic perturbations, formulae in Sec. 5 (Kozai 
1959) can be used for this purpose. 

However, the small divisor disappears in the expres- 
sions of the radius 7 and the argument of latitude L, 
computed by combining the expressions of the semi- 
major axis, the eccentricity, the argument of perigee, 
and the mean anomaly. Furthermore, expressions (13) 
of y and L in the previous paper (Kozai 1959) have been 
proved to hold even for a satellite with very small 
eccentricity. In the present paper itis assumed that the 
eccentricity is a small quantity of the first order. 


ORBITAL ELEMENTS 


Instead of the eccentricity, the argument of perigee, 
and the mean anomaly, the following discussion uses 
the elements \, £, and 7, which are defined as follows: 


A\=M-+e, 
&=e cosw, (1) 
n=—esinw. 


On the assumption that the eccentricity is very small, 
the first-order inequalities of the semimajor axis, the 
inclination, the longitude of the ascending node, and 
the mean argument of latitude, are easily derived: 


da= (A/a) sin’ cos2v, 
di=%(A>/a?*) sin’ cos2n, 
dQ=%}(A»/a?) cosé sin2h, 
dd=}(As/a’)(1—§ sin’) sindd, 


(2) 


where As is the coefficient of the second harmonic of the — | 
earth’s potential as defined by Kozai (1959) on the | 


same assumption; the differential equations for € and y 
are written as follows: 


dé 
—=bn—An sind 
dt 
i] 7 Ay 
— Bul ~~ sin(—\)+- sinan|+C6—, (3) 
4 4 e 
d 


nN il if 
—=—bi—An cosh~ Ba —- Coste ae cosh} 


where 7 is the mean motion and 
b= (A2/a*) (2—3 sin’) n, 
A= (A»/a’)(1—3 sin’2), 
B= (A,/@’) sin’, 
C=3(A3/aAz2) sint. 
The last term in the equation for é is due to the third 
harmonic of the earth’s potential. 
This system of equations is easily solved, if 7, a, and 


n in the right-hand sides are assumed to be constant and 
d is a known linear function of time: 


1 fl 
€ COSw= ey Sinw+ A cosh | = cos(—)) a cossh 
(4) 


€ sinw=é) sinw+A sind 


1 7 
+3: sin(—))-+— sind] +C, 
4 12 
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where e is an integration constant called proper 
eccentricity, and ®=wo+dt. Here wo is another integra- 
_ tion constant. 


COORDINATES 


The radius r and the equation of center »—M for 
this satellite of small eccentricity are written in the 
following way: 


r=al1—e cos(A—a) | 
=al1—e 9 cos(A—w)— A—8B cos2A—C sind ], (5) 
v—-M= 2e sin(A—w) 
= 2e9 sin(A—@)—2B sin2A—2C cosh. 


In the expression (5) a must be regarded as a variable, 
not a constant, and be written as a=G@-+6a. If ap(1— A) 
is adopted as a mean value of the semimajor axis, and 
the third harmonic term is neglected, as was done in the 
previous paper, the inequalities of the radius vector and 
the argument of latitude are written as follows: 


1A, 
ir=a( 1+-— sin cos2h J p 
oe (6) 
6 
Als filial 
aL= -=(--— snt) sin2). 
Dat? 


@ 


It is easily proved that these expressions are identical 
with those of Eq. (13) in the previous paper (Kozai 
1959) and that the definition of @ remains the same on 
the assumption that the eccentricity is very small. 

If the mean values of the radius and the argument 
of latitude are computed with e and w, which are 
derived from the equation 


€ COSw=€) COSH, Ee SiNw=e9 Sinw+C, (7) 


then one can say that instantaneous values of r and L 
are the sums of the inequalities (6) and the mean values. 
‘In other words it has been proved that formulae (13) 
can hold for the small eccentricity. 


MEAN MOTION OF ARGUMENT OF PERIGEE 


If the third harmonic term is neglected, the expression 
of the radius vector is written as 


1A, 
r=a{i—e cos(A—@)-+-- — sini cosh 
nae 


And if the proper eccentricity is very small compared 
with the value of A»/a?, the radius has two maxima and 
two minima. However, even for this case the position of 
the satellite can be expressed by osculating elliptic 
elements. 
In fact, if the new definition of the semimajor axis is 
-adopted, the expressions of ecosw, and esinw are 
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written as 


1 7 

€ COSW= eg const B| cos(—A)+— cossh] 
4 ie (8) 
1 if 

€ Sinw= é9 sina +B - sin(—)) i sinan] + 


The osculating eccentricity is (e sinw)?+(e cosw)?, and 
its constant term is 


er+ (29/72) BP+C?. 
And if the inequality 
eo>2zB+C 


holds, the mean motion of the argument of perigee is 5, 
as in the usual case, and the deviation from the mean 
value @ is written as 


1 7 
€ sin(o=a)=B| sin(—A—@) +— sin(3A—A) 
4 12 


+C cosa. (9) 


However, if the inequality does not hold, and, for 
example, 
(7/12)B> e+ (1/4)B+C, 


the mean motion of the argument of perigee is 37. And 
if C is bigger than the sum of the other three amplitudes, 
the argument of perigee does not have a mean motion 
but oscillates around a certain point. 


DISCUSSION 


At the Smithsonian Astrophysical Observatory, the 
first-order perturbations are taken into account in the 
form of (6) in order to derive the orbital elements of 
artificial satellies. If the eccentricity is small, the 
eccentricity derived is the proper eccentricity and the 
argument of perigee is # when the third harmonic term 
is omitted. However, the proper eccentricity is not 
equal to the mean one for this case, as is easily under- 
stood by the formula (8). When the third harmonic term 
is included, the eccentricity and the argument of perigee 
are represented by (7). 

However, if some disturbing forces other than the 
oblateness of the earth are large enough, one must 
express the variations of the Smithsonian elements as a 
function of time in very complicated way. This condi- 
tion is true for satellite 1960 1.1 (Echo I), because the 
solar radiation effects were very large. 

The Smithsonian eccentricity and argument of perigee 
for this satellite, during the first two weeks after the 
launching, can be expressed quite well by the following 
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expressions: 


e sin = 0.04409 sin (91°99-+2°965#)-+-0.00066 
0.03579 sin (J—2)+0.00188 sin (/-+2 
—0.00519 sin/, 


€ cosw= 0.04409 cos (91°99-+-2°965:) 
+0.03373 cos(i—2)—0.00200 cos(/+2) 
—0.00169 cosi, 


where / is the longitude of the sun, and time is counted 
from August 25.0, 1960 in days. The first terms corre- 
spond to é sine and é cos@, the second term in e sinw 
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Using the Vinti model for an oblate earth it is possible to show that for small values of the oblateness 
parameters the motion of a near earth satellite is almost periodic. To prove this result, bounds on the varia- 
tion of the radial distance and the variation of the angle of inclination are obtained. 


HE purpose of this paper is to extend the results 

of Vinti (1959) on the near earth-satellite prob- 

lem. In a paper that appeared in the Journal of Research 
of the National Bureau of Standards, Vinti showed that it 
is possible to put the force potential into such a form 
that the Hamilton-Jacobi equations are solvable in 
oblate spheroidal coordinates. In the potential form 
for an oblate, axially symmetric earth, we would expect 
not only the even-order Legendre polynomials, but also 
the odd-order ones. Vinti’s potential function provides 
sufficient generality for exact agreement with the first 
three terms of such a gravitational potential. If- only 
even-order Legendre polynomials are permitted, then 
only the first two terms will agree with the Vinti 
potential. Under these assumptions the problem of the 


motion of a near earth satellite perturbed by oblateness’ 


force is reduced to quadrature. 

Throughout this paper, we shall use the notation that 
appears in Vinti’s paper and shall make constant 
references to the various formulae in that paper. 

In oblate spheroidal coordinates, the potential is 


V fala es “14 3(- — R°J»)*r-*"P2,(sind) 
/— — a i ay 2nP, (sin 
[om 2 2 2 


—r+RdJ2 > (—r7R?J. 2)*P. 'on+1(sin§) |, 
6 


where 


C=R J, and J\;=dJz. 
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is due to the odd harmonics of the earth, and the others 
are from the solar radiation pressure. After two weeks, 
the satellite sometimes entered the shadow of the earth, 
so the expressions become much more complicated. 

Of course, one can determine the value of d up to the 
second order by using the present elements (cf. 
Kozai 1955). \ 

The author is grateful to Imre Izsak for his 
suggestions. 
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We will make the assumptions throughout this paper 
that d is between 0 and 10, and we shall also as- 
sume that J2 is between 0 and 10. For the earth 
J2~1.63X10-. 

Vinti found the solution S of the Hamilton-Jacobi 
partial differential equation to be 


Saar f (Pte) 
X [eas + (R+2up+2eip*) (p?-+c7) Pdp 
+f a—9- 


XP —at-+ (b+ 2ndn + 2onc’s?)(1—a2) dp, (57) 


where p=cé and where the upper sign is taken if the 
variable is increasing, and the lower sign is taken if the 
variable is decreasing. 

Because the solution S is separable in oblate sphe- 
roidal coordinates, we can obtain the simple expressions 
for the generalized momenta ,, p,, and pg. 


as _ [eat Et tap aw PtP 


== “(A 
p as Pie ( ) 
= 0S/dq 
—a?+(—k+2 2ac?n?) (1— 
ae P+ (—k+2pdnt 2a?) (1—1?) (64.2) 
iF 
ps=9S/dG=a;. (64.3) 


SATELLITE MOTION 


i If we allow the constants J; and J2 to equal 0, the 
3| coordinate system becomes a spherical one and the 
i) bounds which we obtain from Vinti’s Eq. (57) for the 


SS OA ne = 
= —< = ——— 


two-body problem are 
1/Lo(1+e)<p< 1/Lo(1—e) (65.1) 


r<sinra, (65.2) 


where io is the angle of inclination between the plane of 


the orbit and the polar axis, é is the eccentricity of the 


_ ellipse, and Lp is the reciprocal of the semilatus rectum. 
| We will assume that the initial conditions of the satellite 


are such that R<vp and the satellite would move in a 
Keplerian ellipse if the earth were not oblate. 

We shall see that under the assumption that J2<107 
and d<10~, it is possible to bound the constants of the 
motion in terms of the constants of the two-body 
problem to within “large O” of the oblateness param- 
eters as J. approaches zero (Vinti. 1959). For example, 


if 2 
a1=—p | =r d Abs (—1) "ro 2" 1 R?*J 2" Pon (SinO) 
2ao 1 


dR 
raceme (<r *RI:)*Pans( sin) | ; (66.1) 


rT e 


where do is the semimajor axis of the initial instantane- 
ous ellipse, 79 is the length of the initial radius vector 
and @p is the initial latitude. Hence 


a= —pl (1/2a0)+-O(J2) ], (66.2) 


where 
|O(J2)| <J2/rol2/ (1—J2) ]. (66.3) 
This shows that the total energy of the system 
remains within “large O” of the total energy of the initial 
instantaneous Keplerian ellipse. By “large O” we mean, 
for example, that |a:+(u/2a0)|<CJ2 and for all Js 
sufficiently small where C is some constant. 


Lemma 1. 


From p?=r— (1—7)e and py=r sin@, we can easily 
obtain the following equalities and inequalities: 


0<r—p<e (67.1) 
1>p/P>1—J, (67.2) 
0<9?—sin#<Je (67.3) 
Pt+e=r+er (67.4) 
P+re=r—(1-w)e (67.5) 
(67.6) 


Pe— (r—n)e=P sin’. 


Lemma 2. 


For the initial values of r, 7, p, 6, 7, and 7, 


|ri*| <eoho/(1—eo) 


eoJ 2 
1+ 
: ho 1—€ 
UNS S 
rr 1—Js 
eoJ2 
Se 
2 1-—e) 
| p6| <ho ota eee Ie 
\1—ep r i—J> J 


where h2=p/L. 
Proof. It can easily be shown that 


|*|<eoholo and -|(sin@)*|</o/r?. 
From 777?— (2}—7)C=r* sin’é, it follows that 


mr? — (1— 21) e ]=r7 (sin?6—77) +7? sin8 (sing). 


Hence 
€ohto eoJ 2 
Jo+ho i+ 
es 1—é) in 1—e 
1:7 | Sa a a 
rP—(1-2pP)e fr 1-Ts 
From p2=r— (1—7)e, it follows that 
pp=1T-+ Cm. 
Hence 
‘ €oJ 2) 
| i+ 
' €o 2 1l—e 
| pp| <Ao aa 
—€o r 1—J> 


Lemma 3. 


There exists a constant A depending only on the 
initial conditions such that 


[1+ (k/he)| <AJe (68.1) 
if 

@o<3, AX<22 

eoX<ze, A<10. 


Proof. From Eq. (50) of Vinti’s paper, we obtain 


ee foris < oe 
b= (-+2)( eee 
Pre Pre 
Hence 
eteny\? ’ a Case 
|he+k| = (+e)( : ) ey 
ete ete 


; 


2 2u(p—1)—2( aw'+—+') : 
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It follows immediately that 
C013” K 
= —2u(p—r)— 2(ae't+—') 
p+ 2a 
C? cos*iohe? AT oup? 
prT 2+ 
P+rc r(1—J2) 
I 8 
+—(1-9)e< J oho. 
a 1—€9 
Also 
p?+cn? 2 (1—97)e 2 
a( ) pel1- iG 
pete? P+ 
eo] 2 2 
C 1— 
<he'S 2 pice 
1—@ r 1-J2 
And 
(l-?)epP 
(+e) | — a 
rnc? 
(1-—9)c 2 
= (trini-ctat)| 1] 
rP+7c 
=2(1— 9) (1a)! 
+ (rel 
rnc (7? +-c?)? 
eoJ 2 €oJ 2)? 
i+ 1 
2€0 1—¢ 1—€é 
<JS ahi? , 
ar COM ew 1—J, 
(ee 2 ‘ 
+( ) et). 
1—€) 


Combining the above inequalities, we obtain the 
desired bound. Hence if 
eos a A << 2) 
and if 
eos To A<10. 

For the remainder of the paper, to simplify the results, 
we will assume é9<}. Analogous results are obtained for 
larger values of ép. 

Once we have the bounds on a, and &, it is possible 
to obtain bounds on p and 7 for the full motion. From 
' Eq. (57), we note that both square roots under the 
integrand must remain nonnegative. It is this condition 
which gives us the bounds on p and 7 for the full motion. 
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Theorem 1. 
sin2zy +307 2 
Se 


1+30/2 


To prove this, we note that 
—a3'—R(1—r)-+ (2udn+ 2arcn?) (112) >0 
for all values of 7. Hence 


sin’ip—1? > (1-7?) g(n,J), 
where 
}g(n,J2) | < 30J¢. 
Therefore 


sin’%ip+ (1—7°)30J> 1, 


which implies the desired result. 
From the energy integral, it is easy to show that 


2a | 
Bude eee (71) | 
1— (4Joa0/R) | 


for all time and hence the satellite cannot wander off to | 
infinity. Of course, for this to be true, we must have the | 
condition that R/4ay)>J»2. We can also show that p>R } 
for all time. i 

From Eq. (57) 


O< Cas?+ (p?+c?) (2aip?+2up+k). 


When J/2=0, this equation reduces to 


p 1 
sad Wacker reed bere yea 


a Ly | 
and hence 
esa et (73) 
= a be : 
Lo(1+e0) Lo(1—eo) 
For J2+0, 
2 1 
(e?-+RJ») (-—+2-—) > —S(o,/ 2) 
ao Lo 
where : 
34 
| S(p,J») | <1,| 1+(p?+J2) | szaCt—aaor+—| |, (75) 
0 


where we have chosen the unit of length such that R=1. 
Then if 
Le(1 = e0?) Be ipa 1 
ee (76) 
39— Lo +L? (11+10e9—e?) 


p> 1 for all time. 


[ 
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Theorem 2. 
1—J,1/¢ 1+J,!/¢ 
SS See (77.1) 
Lo(1+eo) Lo(1—e0) 
(1—eo)[ 2eo+J2!/2(1—e) ] 
eS ite (77.2) 
| 67 
where g> 3. 
One might desire to show that 
1—bJ2 1+)dJ> 
ee ps 
Lo(1+eo) Lo(1—e0) 


for sufficiently small values of J2 where 0} is some con- 
stant. This does not appear to be possible. To prove 
theorem 2, we obtain estimates on the amount that the 
roots of Eq. (72) can vary. Since the arguments are 
similar to ones previously presented, the proof will not 


_be included. The value of q is left undetermined so that, 


depending on é», as large a bound as possible on Jz may 
be obtained. 
Once we have shown that p lies between the bounds 


1—Jy}a 14 J,1¢ 
T4160) Palbe.) 


it is then possible to show that the total angular 
momentum per unit mass and the angle of inclination 
are bounded. 


Theorem 3. 
1+e 1—€ 
cos’%y> cos*i> (1—J») cos%%9 (78.1) 
1—J,'a —J pla 
and 


Ly(1+eo) Lo(1—e0) 
> [> —_—__—__ 


1—Jo). 
1—Joi4 eid Wh 1—J a ( 2) 


(78.2) 


as=h cos*tp>=h? cos*7, where h and 7 are, respec- 


_ tively, the total angular momentum per unit mass and 


inclination angle of the instantaneous ellipse. It is 
easily shown that 


ee p/r> 1—Je. 
From the bounds on p 


Lo(1+eo) L(1-e 
—— > L(1+¢ cosy) >————_(1—-J.), 
1—J,ie 1—J,"¢ 
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and hence a bound on L is obtained. From previous 
equations 


L= Lp cos*i/cos%o. 
It then follows that 


(1+e0) ie (1—e9) cos’io 


= => cos*4 > ————————_ ae 2), 
1—J,i¢ 1—J,"ve 


and the desired result then follows. 

The final theorem that we wish to prove is that the 
motion is almost periodic (Goldstein 1956, p. 294). 
Since p and 7 are uniformly bounded, for any value of 
io other than 90°, ~,, p,, and pg are all bounded. 
Furthermore, p and 7 can have extremal values only at 
their absolute maximum and minimum bounds, since 
p and 7 can be zero if and only if 


Leas’ (k+2upt2a1p*) (e+e) ] 
and 


L—as'-+ (—k+2pdn+ 2an0'n?) (1-7?) ] 


are zero. This means that for p and p,, 7 and p,, the 
motion is one of libration and for ¢ and py, the motion 
is periodic. Once this result is obtained, it is obvious 
that the motion of the satellite must be almost periodic. 
However, in general the periods of p, 7, and ¢ will not be 
commensurate. This means that p, for example, will pass 
through its extremal points other than once every 
revolution and that the motion of the satellite will fill 
a volume in space. 

When 7 equals 90°, since the potential is axially 
symmetric, 7 always remains 90°. It is then possible to 
show that the motion is again almost periodic by 
making a change of variables. 

Let sinA=y. The solution of the Hamilton-Jacobi 
equation is separable in the coordinate system p, A, ¢. 
For a polar orbit A varies as ¢ did in a nonpolar orbit. 

pa= (—k+2y6 sinA+2a,c? sin?A)?, (79) 

Hence the motion of A and pa is periodic and this 
implies the almost periodicity of the motion of the 
satellite. 
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The deviations AL and AL’ between the observed and tabular longitudes of the moon and the sun are 
due to two causes. The first cause is the friction of the tides, resulting in a secular retardation of the rotation 
of the earth, and the reaction of this friction resulting in a secular retardation of the true motion of the 
moon. The second cause is the irregular rotation of the surface of the earth due to random currents in the 
interior of the earth. A probabilistic theory of these fluctuations is developed. By a modified method of 
least squares, the coefficients of the secular terms are determined from modern and ancient observations. 


1. INTRODUCTION 


ET AL and AL’ be the corrections to be added to 
the tabular longitudes of the moon (Brown’s tables 
without the empirical sine term) and the sun (New- 
comb’s tables). Both AL and AL’ may be represented as 
sums of secular terms and fluctuations, due to the 
irregular rotation of the surface of the earth. If B* is the 
fluctuation of the moon and T the time in Julian 
centuries from 1900, we may write 


AL=x+yT+2T?+B* (1) 


AL! =x! +y'T+2/T?+ (B*/13.37). (2) 


The factor 13.37 in (2) is the ratio of the mean 
motions of the moon and the sun. For the fluctuations 
the symbol B* is used instead of the usual B, because B 
is defined, according to modern conventions, by 


AL=4"%65+12796T+5122T?-+B, (3) 


and it is desirable to leave open the possibility that the 
coefficients x, y, z in (1) differ from the coefficients xo, 
yo, 20 in (3). In fact, the main problem of this paper is 
the theoretical definition and the empirical estimation 
of the constants x, y, z and x’, y’, 2’ in (1) and (2). 

De Sitter (1927) and Spencer Jones (1939), who pro- 
posed the formulae (1) and (2) and similar formulae for 
Mercury and Venus, originally introduced a factor,Q in 
the B* term in (2), but Jones concluded that the obser- 
vations are compatible with the assumption Q=1. 

Several authors, from Newcomb (1912) to Brouwer 
(1952), have tried to determine the coéfficients in (1) and 
(2) from ancient and modern observations. The method 
was as follows. First, «, y, z were found by minimizing 
the weighted sum of squares of B*. Next, B* was com- 
puted from (1) and substituted into (2), and x’, y’, 2’ 
were estimated by least squares. 

The latter estimation, carried out by Spencer Jones 
(1939) and confirmed by Clemence (1943) led to 
perfectly satisfactory results. The fluctuations of the 
longitudes of Mercury, Venus, and the sun were shown 
to be proportional to those of the moon, and the coeffi- 
cients of T? to be proportional to their mean motions. 
This means that one single clock time correction pro- 
portional to 

2'T?+ (B*/13.37) (4) 


is sufficient to bring the.observed motions of all planets |) 
into-accordance with gravitational theory. a 

The results of the evaluation of x, y, and z were less | 
satisfactory. De Sitter found from ancient and modern | 
observations z=5"22, but Brouwer, giving modern | 
observations more weight, found z=2"22. Both values | 
are incompatible with the theoretical relation, derived 
by Jeffreys (1952); 


z 8.32’ 
—=§.05-- ; (5) 
oi N-+N’ 
in which WV and WN’ are the retarding couples due to the 
lunar and solar tides. On the assumption that tidal 
friction is a linear function of the velocities, Jeffreys 
finds 
N/N'=3.4, hence 2/z’=6.94. (6) 


De Sitter’s value z=5.22, combined with Jones’ esti- 
mate z’=1.23 would give z/z’=4.2, and Brouwer’s 
estimates z= 2.22 and z’=1.01 would give, still worse, 
z/z'/=2.2. Thus, one gets the impression that ancient 
observations contradict modern ones and that both 
contradict gravitational theory. 

Jeffreys’ relation (5) is based on the assumption that 
tidal friction is the only cause of the secular changes in . 
the rate of rotation of the earth. However, in a paper of 
Lord Kelvin (Thomson 1882) it is pointed out that there 
is also a thermodynamic acceleration of the earth’s rota- 
tion, caused by the attraction of the atmosphere by the 
sun. Observation shows that the semidiurnal variation 
of barometric pressure is much greater than the diurnal 
variation. This semidiurnal variation gives rise to a 
gravitational effect analogous to that of the solar tides, 
but opposite in sign. The relation (5) would hold for the 
tidal parts of z and z’ only. As we shall see, the agree- 
ment between theory and observation is slightly im- 
proved by taking the atmospheric effect into account. 

A serious objection to the usual method of determin- 
ing z and 2’ is the principle of minimizing the sum of © 
squares of the B*’s, adopted by all authors. This method 
would be correct if the B* were independent random 
variables, e.g., observational errors. As a matter of fact, 
the B* are highly correlated. Van Woerkom (1953) 
found between values of B in successive years a correla- 
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tion .999, and between values 20 years apart the correla- 
| tion coefficient was still .831. Therefore, another method 
| of estimating «, y, z and their standard errors ought to be 
applied. 

A minor difficulty is the interconnection between the 
estimation of x, y, 2 and of x’, y’, 2’. If we replace one 
estimate of (x,y,z) by another one, we have to change 

the definition of B*, and hence to replace the set 
(x’,y’,2’) by another one. To avoid this, we shall 
introduce instead of (x’,y’,z’) a new set 


p=13.37x'—x, 7=13:312 —3, (7) 


The empirical determination of p, g, and 7 is easy. 
Elimination of B* from (1) and (2) yields 


AL—13.37AL'=—(p+¢T+rT®). (8) 


The left-hand member of (8) is an observable 
quantity. The right-hand member changes very slowly 
with time. Therefore, in observing the quantity on the 
left, the exact clock time is quite unimportant. The 
coefficients p, g, y may be estimated by classical least 
squares, using the full accuracy of modern observations. 
The problem has nothing to do with the irregular 
‘rotation of the earth. It has been solved by Spencer 
Jones and Clemence; their results will be discussed in 
Sec. 3. 

Having determined #, g, and 7, we are left with the 
problem of determining «, y, z, and B* from (1). Here, 
the difficulty is to split the observed AZ into a secular 
part «+yT-+2T? and an irregular part B*. 

This problem will be solved first from the theoretical 
point of view. We shall show that under reasonable 
assumptions the quantity AZ is in fact a sum of a 
quadratic function x+y7T-+27? and a random variable 

_ B* with expectation value zero, which may be repre- 
sented as an integral of another random variable @ with 
respect to time. 

Next, we shall test our hypotheses by means of van 
Woerkom’s correlation analysis. We shall see that just 
such correlation coefficients as van Woerkom has found 
may be expected theoretically. The theoretical con- 
clusion that the differences 6B* are positively correlated 
for small time intervals, but practically independent for 
large time intervals, is confirmed by the data. 

On the basis of this conclusion, a method of estimat- 
ing the unknown parameters «, y, and z will be developed 
and the order of magnitude of the random errors of the 
estimates will be discussed. From this discussion it will 
become clear that there is no contradiction between 
ancient and modern observations, nor between theory 
and observation. Also, there is no need to assume 
changes in the moment of inertia of the earth, as some 
authors have done. - 

The present investigation is mainly theoretical and 
methodical. The numerical values given are only pro- 
visional approximations. To obtain better values, more 
historical investigations concerning ancient observations 
are necessary. 


g=13.37y'—y, 
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Free use will be made of the methods and terminology 
of mathematical statistics. Probable errors will not be 
used, only standard (root-mean-square) errors. The 
square of the standard deviation is called variance. The 
expressions mean value and expectation will be used 
as synonyms. 


2. RETARDATIONS AND FLUCTUATIONS 


We may, theoretically, distinguish between two 
phenomena, viz: 


(1) the secular retardation and the irregular fluctua- 
tions of the rotation of the surface of the earth, caused 
by the friction of the tides and by irregular disturbances 
within the earth; 

(2) the real retardation of the motion of the moon due 
to the reaction of the tides. If the angular momentum of 
the rotation of the earth diminishes, the total momen- 
tum of the earth-moon system remaining constant, the 
angular momentum of the orbital motion of the moon is 
bound to increase, which means that the distance of the 
moon will increase and its angular velocity will decrease. 

Phenomenon (1) gives rise to a correction to the 
apparent longitude of the sun of the form 


AL! =x! +y'T+2eT?+E* (9) 


and to proportional corrections to the longitudes of the 
planets and the moon. The correction for the moon is 


A,L=13.37 (a +y'T+2'T)+13.37E*. (10) 
Putting 13.37E*=B*, we may write instead of (10) 
A, L=13.37 (a! -+y'T+2'T?) + B*. (11) 


Phenomenon (2) gives rise to a correction to the true 
longitude of the moon. If the frictional couple due to the 
lunar tides is assumed to be constant during the 27 
centuries for which observations are now available, this 


correction must have the form 
AoL=—p—qT—rT?. (12) 


The sum of the two corrections is the total correction 
to the tabular longitude of the moon 


AL=A\L+A.L=a+yT+2T?+B*, (13) 
with 
x= 13.37x'—p 
y=13.37y/—q (14) 
2=13.372'—r. 


The formulae (14) show that the quantities #, gq, 1, 
defined theoretically by (12) coincide with those defined 
empirically by (7) in the Introduction. In particular, ris 
the true secular retardation of the moon caused by the 
reaction of the tides. 

It should be possible to determine the retardation r by 
direct observation of the true motion of the moon, using 
an atomic clock independent of the rotation of the earth. 


140 Br 


From observations of this kind, A»Z could be determined 
directly, without observing the sun. Fitting a parabola 
(12) to the observed values of As, the coefficient r may 
be determined. The coefficients » and q depend on the 
adjustment of the atomic clock, but r does not. 

We shall now discuss the determination of #, g, r by 
usual astronomical methods. 


3. EVALUATION OF ~, g, AND 


Spencer Jones (1939) adopted the following set of 
values of x, y, z (expressed in seconds of arc): 


xo=4.65, yo=12.96, 2z=5.22 (15) 
and defined the fluctuation B by (3). Next, putting 
LD’ =4'+y'T+2'T?+ (B/13.37) 


and assuming similar formulae for Mercury and Venus, 
Spencer Jones found 


4’ =1.00; 'y'=2.97, 2/123: (16) 

From (15) and (16) we may compute #, g, and r 
by (7): 

p=8.7, g=26.7, r=11.2. (17) 


Spencer Jones’ remark that z) may be replaced by 
Zot+c with arbitrary c, provided z’ is at the same time 
replaced by 2’+.0748c, is equivalent to the statement 
that r does not depend on the adopted values xo, yo, and 
zo. The same holds for # and gq. 

C. A. Murray, of Greenwich Observatory, had the 
kindness to check the standard errors of z’ and r in the 
evaluation of Spencer Jones. The result of this calcu- 
lation was 


z=1.234.07, hence r=11.2+.9. 


Clemence (1934) calculated z’ anew, making use of 
the 1940 transit of Mercury. The value of 2’ derived 
from this calculation would be z’= 1.07, with a standard 
error of .23. The larger standard error seems to be due 
to the fact that Clemence did not use the older observa- 
tions which Spencer Jones had included. It would be 
highly desirable to make a new calculation based upon 
all available observations. In this paper the value of 
Spencer Jones and Murray, 


r=11.2+.9, (18) 


will be retained. 

If the friction of the tides has not changed since 
antiquity, the value of r cannot have changed either. So 
if z and z’ are determined from ancient observations and 
r computed from (13), we expect to find, within the 
limits of observational errors, the same value of r. 

This expectation is fully confirmed. Fotheringham 
(1921 and 1923), having examined all ancient Greek 
equinox observations, occultations, and eclipses, finally 
arrived at 


2=4.7 and 2=1.52, hence r=15.6. 
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Schoch (1931) obtained a very good representation of . 
all ancient observations by assuming 


z=5.08 and 2’=1.51, hence v=15.00 


(say) .3 in Schoch’s estimate of z’ is assumed, the result- | 
ing standard error of r would be 4 or 5. Hence, there is} 
no contradiction between ancient and modern obser- | 
vations of r. | t 


4, A THEORETICAL MODEL 


Brouwer, at the end of his paper (1952), mentions the | \ } 
possibility of ascribing the irregularities in the earth’s |} 
rate of rotation to “electro-magnetic coupling of the |} 
mantle to a turbulent core.”’ This explanation had been | 
put forward by Munk and Reyelle (1952) and also by | 
Runcorn in Cambridge (England). ‘y 

In the following theory, we need not decide whether | | 
the coupling between the mantle and the core of the }} 
earth is magnetic or mechanical. We start with the i 
following assumptions: 


1. The retarding couple K due to the tides is constant. 

2. The earth’s moment of inertia C is constant. 

3. The mantle is a rigid body. Let ¢ be its moment of | 
inertia and w its angular velocity at time T, hence 


m= Co) 


(20) | 


its angular momentum. Let M be the angular momen- | 
tum of the whole earth, and let the angular velocity 2 | 
of the whole earth be defined by 


M=CQ. (21) 


4. We assume that the friction between mantle and 
core is proportional to the difference Q—w 


Ki=f-@Q-w). (22) 


5. In addition to this friction, we assume another 
couple K», acting on the mantle, due toirregular currents ~ 
in the core. We assume K» to be a random variable, with 
mean value zero, depending on time: 


Ko=V(T). (23) 
The equations of motion are for the whole earth 
Co=—K (24) 
and for the mantle 
co= K+ K —— 
= f-Q—w)+y(T)—K. (25) 
The solution of (24) is 
Q=a—fT, (26) 
with B=K/C. 
To solve (25), we put 
9-6 4(7) (27) 
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daa obtain for ¢ the differential equation 


] 


| cot fo=K'—¥(0), (28) 
ith 
| K’=[(C—0)/C]K. (29) 
i The solution of (28) is 
= y(T—?) 
emo f oval, (30) 
0 c 


vith 


g=f/e and g=K’'/f. 


The solution (30) shows that g=Q—w oscillates 
(bout a mean value go. The deviation is an integral, 
‘lepending on the behavior of the function y in the past 
it times T—i, where ¢ goes from 0 to . If ¢ is large 


compared with 
: to=1/g=<c/f, (31) 


‘he factor e~?‘ is very small. So the important part of 
‘he integral is only that from 0 to 270; the rest may be 
neglected. 

| It is not necessary to assume that the values of the 
random variable ¥(T) for different values of T are 
Independent. It is perfectly possible that a high or low 
value of ¥(T) tends to prevail for a time of the order of 
magnitude 71, say. If this is the case, a low or high value 
of the variable ¢ will tend to prevail for a time of the 
order of 


T=T0+71. 
7 
j 


(32) 


_ Hence, we may suppose that after a lapse of time 27 
the value ¢(7+27) is practically independent of ¢ (7). 
From (26) and (27) we obtain 


\ w=a—BT— 9(T), (33) 


where ¢(T7) varies about a mean ¢go, while values of ¢ 
for widely different times are independent. The mean 
value of w is, according to (33), a linear function of T, 


S=a—8T— ¢o. (34) 


This mean value belongs to a uniformly retarded 
motion. The deviation of 6 from the actual value w is 


0=8—w= 9— Yo 
*y(T—1) 


=— —o'dt, 
0 c 


(35) 


Hence, if our theory is correct, the deviation @ would 
be a random variable with expectation zero. For small /, 
of the order of magnitude r=70+71, we may expect to 
find a positive correlation between @(T) and 6(7+2), 
but for larger ¢ the values 6(T) and 6(T+-2) are expected 
to be practically independent. 

Essentially the same results would apply if c, C, f, 
and K were assumed to be random variables, fluctuating 
about mean values co, Co, fo, and Ko, provided the 
fluctuations are not too large. Even a trend in ¢ or f 
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would not make much difference, provided the rate of 
change of 
K’ G=cK 
i a ae ae 
: agen eae | 
remains small as compared with the retardation 
coefficient in (26): 


B=K/C. 


We shall now compare the theory with the ob- 
servations. 


(36) 


5. THEORETICAL AND OBSERVED FLUCTUATIONS 


From (34) and (35) we find the angular velocity of 
the Greenwich observatory 
Q= wo—BT— é. 


Here wo and 8 are constants, whereas @ is a random 
variable with expectation zero. Integrating with respect 
to T, we obtain the displacement of Greenwich since 


1900 
T 
f wdi=woT —38 =a 
0 0 


The last term represents the deviation of Greenwich 
from a ‘‘mean Greenwich” which has a uniformly re- 
tarded motion defined by the first two terms 7 —38T>. 
Tf this westward deviation 


T 
fea 


0 


T 


édt. (37) 


is divided by the angular velocity #, which for this 
purpose may be regarded as a constant, we obtain the 
clock time correction. If this is multiplied by the mean 
velocity of the moon, we obtain the fluctuation of the 
moon 


xz 
BY=k f adt, (38) 
0 


where & is a known constant. 

To obtain the actual longitude of the moon, we have 
to add to the tabular longitude a correction AZ, 
consisting of 


(i) a correction to the mean motion, of the form 
(Gea 

(ii) a correction AT? for the secular retardation of the 
earth’s rotation proportional to the term 3877 in (37); 

(iii) a correction for the real retardation of the moon, 
caused by the reaction of the tides, of the form (9), 


A.wL=—p—ql—-1T°; 
(iv) the fluctuation B*, given by (38). 


The sum of the four corrections is an expression of the 
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form (1), 


AL=a+yT+2T?+B*. (39) 


We have to compare this with the empirical formula (3) 
defining B, 


AL=xot+yol +2017+ 3B. (40) 
Putting 
x=xota, y=yotb, 2z=2Z0+¢, (41) 
we obtain from (39) and (40) 
B*=B—(a+bT+cT?). (42) 


Hence, if our theory is correct, there exist constants 
a, b, and c such that the B* of (42) is the same as the 
B* of (38). Our problem is now to determine the 
constants a, 6, c, to find out whether the function B* 
defined by (42) behaves like an integral (38), and to 
study the statistical behavior of the random variable 
6(T). 

The determination of a is trivial. Substituting T=0 
in (38) and (42), we conclude 


a=B(0). 


In order to find 6 and c, we form the difference 
quotient 


(43) 


A=6B/8T (44) 


over an arbitrary time interval. If T is the mid-point 
of this interval and 6 the mean of 6 over this interval 
(i.e., the integral, divided by 67), we obtain from 
(38) and (42) 


k= A—(b+2T). (45) 


The mean 6 behaves just like the function 6(T) 
itself: It is a random variable having expectation zero; 
for time intervals 67, and 6T» far apart the means 6(T;) 
6(T.2) are practically independent, but if T: is near Ty 
a positive correlation between 6(71) and 6(T2) may be 
expected. 

From (43) it is clear that the definition of a is a pure 
convention. “Mean Greenwich” was defined so as to 
coincide with actual Greenwich in 1900, that is all. 

On the other hand, 6 and ¢ have a real physical 
meaning: b+-2cT is, according to (45), the expectation 
of an observable difference quotient A=6B/6T. Hence, 
the constants 6 and c may be estimated and the accuracy 
of the estimates judged by the methods of mathematical 
statistics. 


6. THE CORRELATION BETWEEN SUCCESSIVE 
DIFFERENCES 6 B 


If successive observed values of B such as 


Bo=B (To), Bi=B (T1), B,=B(T>), 
were independent, the differences 
SBi=Bi—By, 6Be=B2—Pi, (46) 
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would have a negative correlation coefficient —.707} 
According to our theory, they should have, on thi 
contrary, a positive cotrelation. More precisely, thij 
correlation coefficient should be near 1 if 67; and 67} 
are small compared with r defined by (32), and near 0 i} 
67; and oT» are much larger than Te | 


by theory. Van Woerkom (1953) computed the corel 
tion coefficient between successive 6B’s for one-yea'} 
intervals, and found a serial correlation coefficient i 


7 = .608. 


The true correlation must be still higher, becausy| 
B,—Bo and B,—By-are affected by opposite observa, 


empirical correlation +-.608 strictly excludes the possi} 
bility that the true correlation might be —.707. Hence} 
the method of least squares, applied to the B’s, is not} 
appropriate in dealing with modern observations. 

Van Woerkom also determined the correlation ni 
between one-year differences k years apart, such as 


6B,=B,— Bo and 6Ber1=Beyi— By 
for k=1, 2, --+ up to 20. He found that the correlations] 
diminish regularly from 
71=.608 and 7o=.613 
to 


119= .232 and To9= .040. 


j 


This is in perfect accordance with our theory if we’ 
assume that 7 is of the order of magnitude of 
12 years, i.e., . 
7~0.12. (48) 

Hence, theory and observation combined lead to the 
following conclusion: 
For time intervals 67 of 24 years or longer, successive 
differences 6B are nearly uncorrelated. For smaller time 
intervals, there is a positive correlation between 


successive 6B’s. 


7. DETERMINATION OF @ AND 6 
FROM TWO OBSERVATIONS 


Between 1680 and 1820, the times T for which B is 
available are from 7 to 29 years apart, so we may 
assume the 6B’s to be nearly uncorrelated. From 1820 to 
1952, B is available for every year, but if we pick values 
20 or 30 years apart, we may also treat the 6B’s as if 
they were independent. 

From modern observations only, it is impossible to 
estimate c. Therefore, we shall provisionally adopt 
De Sitter’s value 


2=2)=95.22; hence c=0, 


and find a and 6 from observations made from 1681 to 
1952. We shall see that a and b may be determined from 
the first and last observation only, neglecting all inter- 
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nediate years. We shall estimate the standard error of 

. Next, we shall apply the same method to a larger 

eriod, beginning with the Babylonian eclipse of — 424. 

Ne shall determine a, 6, and c from the best available 

bservations and estimate the standard errors of 6 and c. 
For c=0, formula (45) yields 


A=6B/8T=b+k6. (49) 


The expectation of 6 being zero, the expectation of A 
s 6. An unbiased estimate of this expectation is the 
mean of the observed values of A. Giving the single 
values the weights 67, we obtain the weighted mean 


meer, > OB Br-Bo 
b= = = : (50) 
pei or TT 


fence, the estimate 0 is the slope of the straight line 
irawn from the first to the last point of the graph of B. 

If the 6B are assumed to be uncorrelated and if the 
rariance of 6B is proportional to 67 (both assumptions 
re certainly true for large intervals 67) the weighted 
nean (50) is, according to a theorem of Gauss, the best 
inbiased estimate of 6 depending linearly on the ob- 
erved B’s. “Best”? means having the smallest possible 
rariance. 


8. THE STANDARD ERROR OF 0 


The slope 0, defined by (50), is an unbiased estimate 
f b, i.e., the expectation of 0 is 6. This is seen from (49) 
yy forming the expectation of both sides. 

To judge the accuracy of this estimate, we have to 
letermine the variance of b. Putting = 7,—T, we may 
vrite (50) as 


tb=Bn—Bo=6B\+6Bo+::-+6By. (51) 


The time intervals 67 belonging to the differences 
'B;, 5Bo, --- are all assumed to be of one year. It does 
lot matter that the intermediate values By, ---, Bri 
ire not all known from observation. 

Let o? be the variance of 6B,;. Let p,c? be the 
variance between 6B; and 5B;4,, i.e., the expectation 
wf the product 


RB 267 )(bBs44—0T), 


he expectation of 6B; and 6Bi4, peus béT. From (51) 
ve find the variance of tb 


7 =no?+-2(n—1)pi0?-+2(n—2) 202+ - --+-2pns0”. (52) 


For p;, the true correlation coefficient between 6B; 
ind 6B:;;, we have an estimate 7;., the empirical correla- 
jon coefficient computed by van Woerkom for k=1, 
), +, 20. For k>20 the correlations are, according to 
yan Woerkom, insignificant, but if they were neglected 
iltogether, our estimate of V would become biased. 
Therefore, I have adopted for poi, p22, --- the estimates 


re=e*2 (k= 21, 22, +++). 
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Replacing in (52), the p; by their estimates 7; and the 
factors n—k by n (this is not a bad approximation, 
n being as large as 271), we obtain 


V~no?(1+2r;42re+:--). (53) 


The 7; computed by van Woerkom are the correla- 
tions between the observed differences 6B;=B;—B;_1 
from 1820 to 1950, including observational errors. 
Hence, for o? in (53), we must take the variance of the 
observed 6B,’s from 1820 to 1950. The variance may 
be estimated by 


1 
f= > (6B;—b6T)? (n’=130). (54) 
nN 


Substituting s? for o” in (53) we obtain an estimate of 
the variance of 0: 


ns? 


hare (14+-2r+2ro+- ++). 


(55) 


In (55), 7 is the number of years and ¢ the number of 
centuries of the time interval from which 6 is calculated, 
hence (= .01”. We have 


%=0.35 (56) 
14-2r:+2ret---=22.0 (57) 

Hence we may write (55) as 
= 10/1; (58) 


In our case, the earliest value of B used by Spencer 
Jones is 


Bo=—12.72 (59) 


for 1681, and the latest value computed by Brouwer 
B,=—22.28 (60) 


for 1952.5. These dates yield, by (50) and (58), the 
very inaccurate estimate 


b= —3.5+16.8. (61) 


9. HOW TO OBTAIN BETTER ESTIMATES 


According to (58) the variance of 6 is inversely 
proportional to the time interval. To obtain a more 
accurate estimate we have to take a longer time interval. 

For recent times, the results of measurements by 
atomic time standards are now available. W. Markowitz 
(1959) has compared the cesium measurements with 
the latest moon observations, determining the two 
constants of integration (initial epoch and frequency of 
cesium) so as to provide the best fit with the AZ ob- 
tained by the moon camera for the interval 1955.5 to 
1958.5. He finds a clock time correction At from 
Universal Time to Ephemeris Time (in seconds) : 


1958.0 At=31.10. 
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The time corrections obtained from the cesium 
standard form a smooth function, which can be differ- 
entiated numerically. The derivative of At with respect 
to T (time in centuries) was 58 for 1958.0. 

From At we may calculate B by Brouwer’s formula 
(3) (Brouwer 1958) 


At=24.35+72.316T+29.95T?+1.821B, (62) 


thus obtaining 


1958.0 B=—24.87. (63) 


This value B would yield, by (50), a slightly more 
accurate value of b. 

The accuracy may be increased still more by extra- 
polation. Since Aé and its derivative are known for 
1958, we may find At for 1962 by linear extrapolation. 
Of course, this extrapolation implies a loss in accuracy, 
but this is amply compensated by the longer time 
interval. The result is 


1962 B=—26.041. (64) 


The standard error of the extrapolation is certainly 
less than 1’. To remain on the conservative side, I have 
adopted a standard error of 1’’. 

On the other hand, we may go back beyond 1681, 
Newcomb has computed, mainly from observations of 
Gassendi and Hevelius, the following means: 


1621 B=+24+21 (Eclipses) 
1635 B=— 3+10 (Occultations) 
1639 B=—29+7 (Eclipses) 


Assigning weights 1, 2, and 4 to the three values of B. 


we may form the weighted mean. The result is 
1635 Bo=—14+6. (65) 


The straight line connecting the points (64) and (65) 
has the slope 


(66) 


From three or more modern values of B we could 
determine a, 6, and c by fitting a parabola, but then the 
standard error of the estimate of c would be very large. 
In order to find better estimates, we have to use ancient 
observations. Their larger observational error is made 
good by the much larger factor 7?, by which c is 
multiplied. 


10. ANCIENT OBSERVATIONS 


The following three observations from antiquity seem 
to be more accurate than others: 

1. The cuneiform text CBS 11 901 states that —424 
October 9 (Julian date) a lunar eclipse was seen, 
beginning 10 degrees of time (i.e., 40 minutes) after 
sunset (Kugler 1913, Schoch 1931). Because of the 
short time between sunset and the observation, the 
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time determination is probably more accurate than 
other ancient observations. The sun had set at 17552 
true time Babylon, so the eclipse began, according to tl 
text, at 18532™. Following De Sitter (1927), I assume 
standard error of 10™. I suppose that the observ 
beginning was 3™ after the first geometric contac 
Hence, the first contact was 18529™. The time from tl 
first contact to the middle of the eclipse was, accordi1 
to Ginzel’s “‘Spezieller Kanon” as well as to Schoch 
Neubearbeitung der Syzygientafeln von Oppolzer (Mit 
Astron. Rechen-Inst., Berlin-Dahlem 2, Nr 2), 1250™, 
the middle of the eclipse would be at 20°19™. Accordiz 
to the same tables, the opposition was 2™5 before tl 
middle of the eclipse, i.e., 205165 true time Babylo 
or 17508" mean Greenwich time. At this time tl 
elongation of the moon from the sun would be 180° 
the observation had no error. 

At 5508™ p.m. mean Greenwich time the longitu 
of the sun, according to the table of Newcomb, w 
191°25’30’’. The longitude of the moon, according 
the tables of Brown, was 10°53’20’. Mrs. R. Liist, 
the Max Planck Institute of Munich, had the kindne 
to perform the computation. The difference between t! 
two tabular longitudes is 180°32’10’’. In order to ma 
this difference 180°, the mean longitude L of the mo 
should be increased by 1990’. Hence, we obtain f 
—424 the following equation: 


A(L—L')=1990+305. (6 


From (67), the value of B for —424 may be deduce 
Subtracting (2) from (1), we obtain 


A(L—L’)= (x— x’) + (y—y') T+ (g—2') T?=0.9252B" 


If x,y, 2 and x’, y’, 2’ are replaced by Spencer Jon 
estimates, we have to replace B* by B. Hence, B may 
computed for —424. The result is 


Year —424 B=70+330. (6 


2. Ptolemy used (Almagest IV 11) a lunar eclipse 
—382 December 23. The beginning of the eclipse w 
observed at Babylon, according to Hipparchus quot 
by Ptolemy, half an hour before sunrise. Oppolzer ai 
other authors supposed that the eclipse had be 
observed in Athens or Ionia, and that Hipparchus h. 
reduced the observed time to Babylon, but I have shoy 
elsewhere (van der Waerden 1958) that this is impc 
sible. The time indication “half an hour before the e1 
of the night” may not be very accurate, but in any ca 
the eclipse must have been visible before sunrise 
Babylon. 

Sunrise was 659™ true time Babylon. By “half ; 
hour” Hipparchus probably meant half an hour 
night, or 35™, but he may also have had in mind 30 
The original Babylonian record might have stated th 
the eclipse began 7 or 8 or 9 degrees of time (i.e., 28 
32 or 36 minutes) before sunrise. Therefore, I sh: 
adopt 32 minutes before sunrise, or 6"27™ as the time 


_ observation, and 6"24™ true time Babylon as the time 
| of the first contact. If we had the Babylonian record 
| before us, we could assume, as in the case of the eclipse 

of —424, a standard error of 10™, but since we have only 
a third-hand record, we shall adopt a standard error of 
_ 12. Thus we obtain the value 


Year —382 B=+1010+400. (69) 


The difference between the values (68) and (69) is 
less than twice its standard error. 

3. According to Schoch, the most important ancient 
observation is the contact between Spica and the moon 
observed by Timocharis —282 November 8/9 in the 
middle of the 10th hour of the night, as the moon had 
just risen from the horizon (Almagest VII 3). Schoch 
(1931) assumed that the occultation took place half an 
hour after moonrise. From this assumption he derived 
the following correction to the tabular longitude of the 
moon according to Brown’s tables: 


AL=—3.444+2.70(T+1)+5.08(T-+1)? 
=4,344-12.867+5.087°. (70) 


By this correction, he obtained a very good representa- 
_ tion of all ancient occultations and eclipses. Neugebauer 

(1931) showed that also a Babylonian occultation ob- 

servation of —418 is well represented by Schoch’s 
~ elements. 

For —282, Schoch’s correction (70) would give 
AL=+2130”. By (3) this implies B=—65”. If this 
value were derived from the occultation of — 282 only, 
its standard error would be something like 16™, or 530’; 
but since Schoch’s elements also agree with other ob- 
servations, we may diminish the estimated standard 
error to 12™, or 400’’. Hence, we obtain 


Year —282 B=—60+400. (71) 


Fotheringham and De Sitter attached great import- 
ance to the solar eclipse of —128. It was believed that 
this eclipse had been observed with great accuracy by 
Hipparchus. However, the reports of Cleomedes, 
-Pappus, and Theon do not say that Hipparchus has 
observed this eclipse, but only that he used a certain 
eclipse in his calculations. Neugebauer (1929) has 
pointed out that the eclipse Hipparchus used may well 
be the one of —309. So we have to discard the “eclipse 
of Hipparchus.” 

Combining the values (68), (69), and (71) with 
weights 3, 2, and 1, we obtain a weighted mean 


Year —386 B=+360+220. (72) 


Sachs’ publication of late Babylonian texts (1955) 
contains many more eclipse records from —730 to —64 
(at least). Probably a better estimate of B for —400 may 
be derived from these records. Also the “lunar sixes,” 
i.e., the observed time differences between rising and 
setting of sun and moon, may be used for the same 


- purpose. 
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11. DETERMINATION OF a, 6, AND ¢ 


From now on, we shall drop the notational distinction 
between the true values a, b, c and their estimates G, b, 
and ¢. We shall denote any estimate of c by the same 
symbol c, as usual in astronomy. In order to estimate the 
constants a, b, and c, we have to combine the ancient 
observations with medieval and modern ones. 

Ibn Yunis has collected 25 carefully recorded eclipse 
observations made between 829 and 1004. The observa- 
tions have been reduced by Newcomb (1879) and Nevill 
(1906). From Newcomb’s results Brouwer (1952) com- 
puted a clock time correction for 950 


F=1.821B=—511*+ 2308. 
This gives us a new point of the fluctuation curve 
Year +950 B=—280"+120”. 


Our problem is now to fit a parabola to the four 
points 


Year —386 B;=+360+220 
+950 B.=—280+126 
+1635 Bi=—14+6 
+1962 Bo=—26+1. 


The usual least-squares method to fit a parabola 


B=a+0bT+c?? (73) 


must be modified in our case, because the variables B; 
are dependent. According to (42) we have 


B=a+6T+cT?+ B*, (74) 


where the B* are random variables with expectation 
zero. We shall adopt the convention 


Byo*=0 for 1962. 
If there were no observational errors, the differences 
ui= Bs —By_i* 


would be nearly independent random variables and the 
B,* would be given by the formulae 


Bot=0 

B= 
B*=uU,+u2 
B= Uy+ Us U3. 


To these sums we have to add the observational 
€ITOTS Vo, V1, V2, 13 Adding also the quadratic terms in 
(74), we obtain 


Bo=a+bTot+cT i+ 
By=a+bT14+-¢cT Yt m+ 
Bo=a+bT2+¢T P+ m+ 02 
Bs=a+bT3tcT + ust ust ust. 


(75) 
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The wz; and 2; are independent random variables with 
expectation zero. The variance of #; is, according to (58), 


s2=770(T;1—T;). (76) 

The variances of v9, v1, v2, vs are 

oe=12, of=6, of=126, o°=2207, (77) 

Putting X;=B;—B;1, we find 

X1=6(1T,—Tp) te(T2—T?) +m+01— 00 
Xo=b(T2—T)) te(T2—TY)+usti2—2; (78) 
X3=6(T3— Te) +e(T3—T 2?) +us+03— 02. 

The variances and covariances of the X; are 
Qy=SPtortor ay=—acr 
da=SY¥torvtor a3=0 (79) 
Q33=S3+ost+ox da=—o2. 


In our case @12 is so small as compared with ay; and 
2, that it may be replaced by zero, which simplifies the 
calculation. Besides, #; and #2 are not strictly inde- 
pendent but have a small positive correlation, so that 
putting aj2=0 is even an improvement. The problem is 
to find best unbiased estimates of 6 and c from the 
observed values X1, Xo, X3. To simplify this problem, 


we shall replace the X; by uncorrelated variables 
Yi=Xi 
Y2=aX1+X2 
V3=BXityXot Xz. 


In our case, we may choose a=8=0. The condition 
for y is that Y2 and Y3 are uncorrelated, i.e., that their 
covariance is zero. This gives 


Ya22-+ d23= 0, 


hence 
Qf — 03/A22= 0.75 (80) 
and 
VY, = Xy 
Y2=X2 (81) 


Woes 0.75X0+- X3. 


To the uncorrelated variables Y; we may apply the 
usual methods of least squares. The observational 
equations are obtained from (78) by substituting the 
observed values X; on the left and omitting the terms u; 
and 2; on the right. The independent observational 
equations resulting from the transformation (81) are 


—3.27b6+6.64c=VYi;=+12 (weight 24.4) 

—6.856+83.2c=V2=—266 (weight 3) 

—18.50+ 495c=Y3=+440 (weight 1). 
The normal equations are 


744b—11 397c= —3631 
— 11 3976+-266 870c= +153 350. 
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The solution is 


b= 11.3415.6 


82) | 

c= 1.06+0.82. ( | 
The coefficient a is determined by the condition } 
Bo*=0 for 1962, or | 


Bo=atbTotcT? (To=0.62). 


This gives | 
t= 2970) (83) | 


12, DETERMINATION OF (#,¥,2) AND (2',y',2’) 
Having found a, 6, and c, we might compute x, y and | 
z from (41): 


y=yotb, z=z0+e. 


This would be the right procedure if we had deter- | 
mined a, b, and c from observations of the moon’s | 
longitude L. However, most ancient and medieval | 
observations used in the preceding section give the | 
moon’s elongation D= L—L’ only, and the sun’s longi- | 
tude is not known accurately for ancient times. There- 
fore, we have to adopt a slightly different procedure. 

The deviation of D from its tabular value is | 


AD=A(L-L')= (a—2#')+ (y—y)T 


x=xXo+a, 


+ (g—2!)T?-+0.9252B*, (84) 
Putting | 
x—x’=u, y—y'=2, 2-2/=v, (85) 
we may write for (84) 
AD=u-+0T+wT?+0.9252B*. (86) 


All ancient and modern observations may be con- 
sidered as observations of AD. We may determine wu, 
v, w from these observations by the method of the 
preceding sections. Starting with provisional values 
Uo, Yo, Wo found by subtracting (16) from (15), 


Uo= 3.65, %=9.99, wo=3.99, 
we may define B by 

AD= up tv0F +woT?+0.9252B 

and determine w, v, w by fitting a parabola 

B=a+6T-+cT? (87) — 


to the observed values of B. The final estimates of w, 2, 

w are now found from 
u=uo+0.9252a 
V= Vp + 0.92526 
w=wot0.9252c. 


(88) 
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Substituting (82) and (83) into (88), we find 
u= — 26.0 
v= 20.4+14.4 
w=4.97+0.76. 


The standard errors of v and w are found, of course, 
by multiplying the standard errors of b and c by 0.9252. 

We now combine (89) with (17). Solving the 
equations 


(89) 


x—x’=u 13.37x'—x=p 


y—y'=0 13.37y'—y=q (90) 
2—2=w 13.372'—z=r, 
we obtain } 
=—274 x’ =—1.40 
y= 24.2 yy’ =3.81 (91) 
2=6.28+.82 2’=1.31+.10. 


The value of z’ does not differ much from Fothering- 
ham’s estimate 2’=1.52 derived from Greek equinox 
observations and solar eclipses. 


13. COMPARISON WITH GRAVITATIONAL THEORY 


Jeffreys’ theoretical relation (5) holds for the tidal 


parts of z and 2’ only. Let us put 
Z= 21+ 22 
3! =21' +20! a 


where z; and 2,’ are due to the tides, z2 and 2’ to the 
atmosphere. 

Lord Kelvin’s estimate of the time correction due to 
the atmosphere was 28777, but a more accurate com- 
putation based upon Simpson’s data (see Holmberg 
1952) gave 2:947°. The motions of the moon and the sun 
during this time are 


otre and -0.127°. 


Hence, we have 
z=—1.61 and 2/=—0.12. 
The estimates (91) for z and 2’ now yield 
21=2—22.= 7.89+0.84 
ay =2/—29'=1.43+0.11. 


The quotient is 
Bist = 9.5. (93) 


According to Jeffreys’ theory, the ratio 2:/z:’ ought 
to lie between 6.4 and 6.94, probably nearer to the 
latter value. A reasonable estimate seems to be the 
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weighted mean of 6.4 and 6.94 with weights 1 and 2: 
21/21 =6.76+0.10. (94) 


Within the limits of observational accuracy, (93) 
coincides with (94). So there is no contradiction 
between observation and theory. 

The theoretical relation (94) may be combined with 
the observational relations 


2—2/=w=4.97+0.76 
13.372’—z=w’ =11.2+0.91 


(95) 
(96) 


in order to find, by the method of least squares, better 
estimates of z and 2’. Still more accurate estimates may 
be obtained by using Babylonian eclipse records (see 
end of Sec. 10), recent Mercury transits (Sec. 3), and 
atomic clock observations (end of Séc. 2). 

Recently, Holmberg (1952) put forward the hypoth- 
esis that there is no retardation of the earth’s rotation, 
the acceleration caused by the atmosphere being nearly 
equal to the retardation caused by the friction of the 
tides. In this case, z’ would be zero. By (7), this would 
imply z=—r, hence z would be negative, and w=z—2' 
would be negative too. This is definitely excluded by 
our empirical results 


w=4.97+0.76 
z= 1.31+0.10. 
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The analysis of the joint distribution of distance, redshift, and apparent magnitude of field galaxies shows 
that, under certain conditions, the dispersion o of redshift of field galaxies can be estimated empirically. 
This is possible by an indirect method, through the study of regression of certain functions of the redshift 
on the apparent magnitude. For example, the regression of the square of the redshift may be used for 


this purpose. The method has been applied to data published by Humason, Mayall, and Sandage. The 


estimate of the variance o? obtained is “negative zero,’ namely — (66 km/sec), with an estimated standard 
error of (121 km/sec). This result seems to indicate that the dispersion of redshift of field galaxies is small 
and, thus, the redshift of field galaxies is a reliable distance indicator. 


INTRODUCTION 


HE redshift of galaxies belonging to clusters must 
show dispersion whether or not the Hubble law 
of redshifts reflects the expansion of the universe. For 
the durability of a cluster, its members must be in 
relative motion and the Doppler effect of their peculiar 
velocities will produce an observable dispersion of 
redshift. The empirical estimation of this dispersion is 
connected with the notorious difficulty of distinguishing 
between true cluster members and background or 
foreground objects. Recent estimates of the dispersion 
(Oort 1958; Burbidge and Burbidge 1959) are of the 
order of 1000 km/sec. 

The question of the existence of dispersion in redshift 
of field galaxies, and of the order of magnitude of this 
dispersion if it does exist, is of considerable interest from 
several points of view. Depending on the numerical 
value of this dispersion, the redshift is or is not a good 
distance indicator. Similarly, the probability of chance 
collisions between field galaxies (Page 1961) depends on 
the magnitude and frequency of deviations in their 
velocities from a strict realization of the Hubble law. 
Finally (Ambartsumian 1959) any attempt at the 
construction of a comprehensive stochastic model of the 
realm of galaxies must involve a more or less. specific 
assumption regarding the redshift of field galaxies: The 
redshift is a single-valued function of distance or is a 
random variable with conditional distribution depend- 
ing on distance, etc. 

Twenty-five years ago Hubble (1936) recognized the 
importance of the problem of dispersion in redshift of 
field galaxies and produced the first estimate. So far as 
we know, Hubble’s estimate of the dispersion, amount- 
ing to about 200 km/sec, remains the only estimate of 
this parameter. It was arrived at by plotting the redshift 
Z of 29 presumed field galaxies against their estimated 
distances £ Of the 29 points of the scatter diagram, 
28 formed a more or less convincing group along a 
straight line through the origin. The twenty-ninth 
point was an outlier which Hubble decided to neglect. 
The mean square of deviations of the 28 retained 
points from the fitted regression line passing through the 
origin was taken by Hubble as representing the total of 


the_several sources of variability, namely the intrinsic } 
dispersion of the redshift of field galaxies at fixed | 
distance, the effects of errors in estimating the distances | 
of particular galaxies and the errors in measuring | 
redshift. While the total dispersion estimated from the } 
mean-square residual was 260 km/sec, Hubble thought } 
that the intrinsic dispersion of the redshift must be of | 
the order of 200 km/sec. | 

As is well known, since the time of publication of | 
Hubble’s paper the confidence in the early estimates of | 
distances of individual galaxies declined considerably | 
and the general distance scale has been repeatedly | 
revised. This circumstance poses the problem of devising | 
a method of estimating the dispersion in redshift based | 
on only such variables related to galaxies that are open | 
to direct measurement. The redshift Z and the apparent | 
magnitude m of field galaxies are examples of such © 
variables. Also, one might think of properly defined — 
visible diameters. 

The purpose of the present paper is to investigate the 
possibility of estimating the conditional variance of the 
redshift Z of field galaxies, for fixed distance é, using 
the measurements of Z and of the apparent magnitude 
m. For obvious reasons, such estimation, if at all 
possible, must be indirect. The basic assumption under- 
lying the theory given below is that the stochastic — 
process representing the universe is stationary not only © 
in the three spatial coordinates but also in the fourth ~ 
time coordinate. In other words, it is assumed that (i) 
either the universe is static, with the redshift represent- 
ing the effect of a factor other than the velocity, or 
(ii) that the universe is expanding but in a steady state 
or (iii) that the available data on redshift and magni- 
tude refer to bright galaxies, located at relatively small 
distances, so that the effect of temporal nonstationarity 
of the stochastic process considered may be neglected. 
The other basic assumptions are those given in the 
earlier publications (Neyman 1955; Neyman and Scott 
1952). In particular, we are limiting ourselves to con- 
sideration of Euclidean space. Further assumptions are 
enumerated below. 

While the present paper is limited to the consideration 
of a stochastic process stationary in all four dimensions, 
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we wish to emphasize that the hypothesis of stationarity 
in time is no more than exploratory and that it is hoped 
to consider the alternative hypothesis also. In fact, 
\preliminary results obtained in this direction indicate 
‘that there are additional difficulties in estimating 
jthe dispersion of the redshift under the alternative 
‘hypothesis. 

| On the assumptions just indicated, if there is no 
dispersion in the redshift Z, then the equation 


Z=Hé ay) 


bids, where H is the Hubble constant. If there is a 
dispersion, then it is natural to assume that the condi- 
‘tional expectation of Z given é is 


E(Z|§)=HE. (2) 

The conditional variance of Z given £ will be denoted by 
o°(&)= EL (Z—HE)| €] 

= E(Z?| £)— Be, (3) 

Our problem is to estimate o?(£) empirically. In order to 

do so, admitting the possibility that o7(&)>0, it is 

“necessary to make some assumption regarding the 


“nature of the functional relationship between o°(é) and 
_£. The natural thing is to write a tentative series 


| 0%()=o-+ab-+BE+-: 


(4) 


where a”, a, 8, --+ are undetermined coefficients, and 
then to try to estimate them in turn. As will be seen 
below, in principle this is possible but practical diff- 
lalties are considerable. For this reason, we limit our- 
selves to a closer study of two possibilities, One is that 
_o°(€)=o?=constant for all ~ The other possibility, 
which on intuitive grounds seems plausible, is that the 


_ dispersion of Z is proportional to distance é so that the 
variance is 


PE)=RE, (3) 


where / is a Hubble-like coefficient of proportionality. 
Here, then, the estimation of o?(£) reduces to the estima- 
tion of he. 

Our results include the following oes points: 

(i) On the assumption that the conditional variance 
of Z given & is proportional to the square of £ as in (5), 
it is found that the estimation of h is possible only if one 
knows exactly, or if one assumes that one knows, the 


luminosity function of field galaxies of at least one 


morphological type. As things stand now, reliable 
attempts at the determination of the luminosity func- 


_ tion are based on treating Z as a distance indicator and 
on using (1). Thus it appears that, if one adopts the 
hypothesis (5) with the possibility of having /?>0, 
attempts at an empirical estimation of / must involve 


a vicious circularity: First we assume that o(&) is 
zero and use (1) to estimate the luminosity function of 
field galaxies; next we use the luminosity function so 
obtained in order to estimate # on the assumption that 
o(£) need not be zero. 
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(ii) On the assumption that the conditional variance 
of Z given é is independent of &, say 


a°(£)=0", (6) 


the empirical estimation of o? is possible, without mak- 
ing any specific assumptions regarding the luminosity 
function of the galaxies concerned and without circu- 
larities. The prerequisites for such estimation are 

(a) The availability of a reliable catalogue of field 
galaxies (that is of galaxies forming ‘‘clusters” of one 
member) compiled with an effort at some sort of 
completeness conditioned by the instruments available. 
This catalogue must give measurements of the redshift 
Z and of the apparent magnitude m of each galaxy. 

(b) Another important assumption is that there be 
available a morphological classification of field galaxies 
fine enough to ensure that for any two galaxies classified 
as belonging to the same type the relative ease of 
measuring Z and m depends upon the apparent magni- 
tude of these galaxies and on nothing else. To illustrate 
this point: An irregular galaxy Gi and a concentrated 
elliptical G,; may well have exactly the same apparent 
magnitude, but the ease of measuring their redshift may 
be very different. 

(c) The calculations are based on relatively bright 
galaxies so that the effect of dimming due to redshift on 
the ease of measurements may be neglected. To illus- 
trate this point: The corrections of magnitudes for the 
effects of redshift given for field galaxies in the recent 
HMS catalogue (Humason, Mayall, and Sandage 1956) 
almost never exceed 0.12 mag. The hypothesis made 
amounts to the assumption that within the range of 
magnitude covered by this catalogue the dimming of 
the order about one tenth of the magnitude does not 
influence materially the decision as to whether to 
include a given galaxy in the observational program 
or not. 

(iii) The theory developed is applied to measure- 
ments of redshift and apparent magnitudes of field 
galaxies as published by Humason, Mayall, and 
Sandage. Because of the similarity of luminosity func- 
tions, the galaxies included in the study were divided 
into three categories: faint, medium, and bright. For 
medium-bright field galaxies the estimate of dispersion 
of redshift, including imperfections of measurements 
and supposed the same at all distances, is 121 km/sec. 
For the other two types, the faint and the bright, the 
estimates of the variance of redshift come out negative. 

The estimates of o? obtained are understandably 
shaky but it is hoped that, with the increase of observa- 
tional data and, particularly, with the development of 
finer systems of morphological types (van den Bergh 
1959, 1960; de Vaucouleurs 1959) more reliable esti- 
mates may be obtained. 

The results obtained are connected with the general 
theory of the spatial distribution of galaxies developed 
earlier (Neyman and Scott 1958; Neyman, Scott, and 
Shane 1956). 
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NOTATION AND SPECIFIC ASSUMPTIONS 


Our basic assumption is that we are dealing exclu- 
sively with field galaxies, that is galaxies Poisson-wise 
distributed in space. For each such galaxy we shall 
consider its redshift Z, its apparent magnitude corrected 
for redshift m and its distance &. 

Another assumption is that there is developed a fine 
morphological classification of galaxies with the follow- 
ing property. For each type of galaxies, symbolized by 
the letter /, there exists a probability $:(m), depending 
only upon the apparent magnitude m that a galaxy of 
type ¢ will be included in the catalogue. Thus, any two 
galaxies of the same type located in the part of the sky 
covered by the catalogue and having the same apparent 
magnitude m have the same chance of being included in 
the catalogue. However, for two galaxies of different 
types ¢’ and #’”’, respectively, the corresponding proba- 
bilities @(m) and d-(m) may well be different. No 
other hypotheses are made about the probabilities 
¢:(m). The probabilities ¢:(m) will be called selection 
probabilities. 

All the theory developed below refers to galaxies of 
just one isolated type. For this reason further references 
to types and the subscript ¢ will be omitted. They will 
reappear in the section giving the empirical results. 

For the type of galaxies considered we use the symbol 
pm (MZ) to represent the appropriate luminosity func- 
tion, that is, the probability density of the absolute 
magnitude. Considered as a random variable, the 
absolute magnitude is denoted by 9. They symbol M 
stands for a particular value of 9. Considered as 
random variables, the apparent magnitude of a galaxy 
and its distance will be denoted by uv and 3, respectively. 
Their particular values will be denoted by m and é. 
We have 

w=I—5+5 logicg (7) 


or, denoting a=5 logue 
w=I—5-+a loge. MAS) 


Without making any restrictive assumptions, the 
distance © has a degenerate distribution. However, we 
shall be concerned with conditional distributions, given 
that a galaxy has been selected for inclusion in the 
catalogue. The probability densities deduced under this 
restriction will be marked with asterisks and described 
as the “catalogue densities’ of the variables considered. 

It is known that (Neyman and Scott 1959), whatever 
be the luminosity function and whatever be the 
selection probability for a given type of field galaxies, 
the joint catalogue density of the apparent magnitude 
and of the distance is represented by the formula 


Pu,z*(m,£)=Co(m)Ppm(m+S—alogt), (9) 


where C is the norming constant. Formula (9) is basic 
for all that follows. 
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CATALOGUE REGRESSIONS OF CERTAIN FUNCTIONS 
OF REDSHIFT ON THE APPARENT MAGNITUDE 


The whole of the present study is based on an aspect! 
of the magnitude redshift relation. We begin with 
formula (2) giving the conditional expectation of Z: 
given £ and supplement it by 


E(2| &)=HP+0°(é). 


Also, if one supplements the hypotheses made about 
the distribution of redshift by an extra assumption 
that, for a given &, the distribution of Z is normal, one 
can write 


(10) 


4 


E(Z3| £)= H+ 3H £02(£) (11) 


E(Z4| £)= H+ 6H #0%(E)+304(). (12) 


Since ~ is unobservable, the chances of estimating’ 
o°(é) depend upon the conditional expectations of 
functions of Z given u=m. These can be obtained by 
using (9) and deducing the conditional distribution of 
= given u=m. We have 


(13) 


pz*(é|m)= Cif pn (m+S5—a loge) 
with 
cris f 2p (m+5—a logé)dé. (14) 
0 
According to a familiar formula, 
E(Z*|m)= E{ELZ*|2 || m}. (15) 


Also, using (13), 
ELZ*| m]=Cy ‘f +p (m+5—a logt)dé. (16) 
0 


Introducing a new variable of integration 


t=m+5—a logé (17) 
and factoring out terms independent of ¢, we find 
EEF | m]=CreGt® mtorlag,,, (18) 
where 
+00 
Qr= f e— G+h) tla hoy (d)dt (19) 


—3O 


is a factor depending upon the luminosity function of 
the galaxies considered. Noticing that E(=°|m) must be 
equal to unity, we obtain the following formula: 


E(E*| m) = (gx/ poet !4= (pe/ po) LOR mt 8)/8, 
Using (2), (10), and (15) this leads to the expressions 

E(Z|m)=H(g:/ e105 410" (say) (21) 

E(2?|m)= H?(¢2/ go) 10? +5)/5+ El o2(&)|m]. (22) 


If one assumes that o?(£) is proportional to #, then 


(20) 
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(22) reduces to, say, 
Ey(2#|m)= (PHI) (go/ e128, (23) 


Alternatively if one assumes o”(§)=0?=constant, then 
formula (22) gives, say 


E,(Z?| m) = H?(g2/ po) 10?" )/5+ 9? (24) 
= Bi0?/5+¢2 (say). (24) 


Formulae (20), (23) and (24) appear of considerable 
interest. They may be summed up in the following 
proposition. 

Basic proposition. li the classification of galaxies is 
fine enough to ensure that the probability of a galaxy of 
a specified type being included in the catalogue depends 
upon this galaxy’s apparent magnitude (corrected for 
redshift) and on nothing else, then for the catalogue 
galaxies of a single type: 

(a) The regression (21) of redshift on (say) «= 10/5 
must be linear with intercept zero, with a slope A 
depending upon the luminosity function. 

(b) If the conditional variance o?(£)=/?2, then the 
regression (23) of Z? on x? is also linear with intercept 
zero. Its slope depends both on /? and on the luminosity 
function of the galaxies considered. Thus, as long as 
the luminosity function is uncertain and no extra 
assumptions are adopted, no empirical study of magni- 
tude and redshift of field galaxies can distinguish 
whether 4>0 or h=0, 

(c) If the conditional variance o?(£)=o?=constant, 
then the regression of Z? on x, given by (24) is linear 
with a slope depending upon the luminosity function of 
the galaxies and with an intercept exactly equal to o”. 

It follows from (c) that, without making any assump- 
tions regarding the luminosity function, the estimation 
of o” is possible and reduces to the estimation of the 
intercept of a linear regression with an uncertain slope. 

Actually there is a slight advantage in considering the 
regression of 


Y=(Z—Ax)?, (25) 


where A is the least-square estimate of-A. The reason is 
that the regression of Y is of a form similar to that of Z? 
but the conditional variance of Y is somewhat smaller 
than that of Z*. The conditional variance of Z?, com- 
puted on the assumption of normality, is 


= E(Z4|m)—E*(Z?|m) 


= 10H [ (4/90) — (2/0) fe! 
+400? (¢2/¢o)o7x?+204. (26) 


In order to compute the variance of Y, we must first 
decide on the weights to be used in computing A. 
Because of the suspicion that o? is not a very large 
number, and because the property of unbiasedness of a 
least- ~square estimate is not affected by using wrong 
weights, i in computing A we act as if o? were zero and 
use as weights the values of a. 

Suppose we have » field galaxies of a specified 
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morphological type and let «;=10”*/* and Z; refer to the 
ith of them. Then easy calculations give 


1 
4=- > Zi/x; (27) 
Wt 


and the unbiased estimate of the variance of A is, say 
6={)0 Z2/x2—-A DL Z;/x3/n(n—1). (28) 
Further easy calculations give 


n—2 


E(Y;|m,) =o? 


n—2 1 n 
+x2] 3 ape Lio +Bsa)/xa | (29) 
n neil 


In other words, the regression on x? of Y; defined by 
(25) with (27) is again linear with slope depending on 
B and, therefore, on the luminosity function, and with 
intercept practically equal to o?. Neglecting the terms 
of the order of 2+ and higher, and also the terms 
proportional to powers of o”, the conditional variances 
Z2 and of Y;, given m;, are easily found to satisfy the 
relation 


Var(Y;)= teen —— (= 2) 
gies igo 1\? 
-ef2224(2)] 
Pol Yo Yo Yo Yo 


=Var(Z2)—4E(E| m:) 


X E{Z(2—E(Z|m,) P| m;}. 


The last line of formula (30) indicates that the con- 
ditional variance of Y is necessarily less than that of 
Z2. The first line shows that in the efforts of estimating 
empirically the first term of formula (29) the particular 
values of the Y; should be assigned weights Q;=4;—*. In 
applying the least-squares procedure we notice that, 
because the estimate A intervenes in all the Y,, the 
latter are not strictly independent. However, if the 
number z of observations is substantial, the covariances 
between these variables are small. Neglecting them, the 
least-square estimate of o? is found to be, say 


(30) 


nN is a3 
oF ba}, (31) 


n— 


where b stands for the least-square estimate of the 
coefficient of x2 in the right side of (29), namely, 


pa (x2—-2)0.Y; 


= (32) 

ys (x2— x)20; 
w= SE ey CINE Fai (33) 
Var tino, i+. (34) 
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The leading term in the estimate of variance of (30), 
based on the assumption that the true value of o? is zero, 


is, say 
on) ea 
n—2 Sey 
+ (35) 


If o? is estimated separately using several different 
types of galaxies and if ¢? and S? stand for values of 
(31) and (35) computed for the fth type, and if one 
makes the plausible assumption that o? does not depend 
upon the type of field galaxy, then the several inde- 
pendent estimates of o? can be combined into one, say 


ap Gi Spaces (36) 
The asymptotic variance of this estimate is, say 


SP= (Le S74. 


EMPIRICAL RESULTS 


The above formulae were used to obtain estimates of 
the conditional variance of redshift of field galaxies for a 
fixed distance. The material used was that in the HMS 
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Fic. 1. Scatter diagram and | 
fitted regression line of Y on m. 
Faint galaxies. 


catalogue. In order to comply as far as possible with the 
requirements regarding the postulated mechanism of 
selection, several galaxies were eliminated which, ac- 
cording to the information given either in the paper 
quoted or obtained directly from N. U. Mayall, were 
selected for some special reason. Also an effort was made 
to eliminate known cluster members. 

Because of the dependence of the relevant regressions 
on the luminosity function of the various types of 
galaxies, the earlier theory (Neyman and Scott 1959) 
was used to obtain estimates of the several luminosity 
functions, This was done for types of the original 
classification published in the catalogue. In each case it 
appeared that the normal frequency distribution is not 
an unreasonable approximation to the distribution ob- 


Taste I. Constants of selection probabilities and 
of luminosity functions. 


Group Type a B Mo om n 
Faint Sc, Scp, SBc 11.6 0.75 =15.4 1:28 94 
S0, SOp 11.7 1.10 | —15../) Sie25eanoo 
E4-E7, Ep 10.9 1.30 "=15..7 esses 
Medium Sb, Sbe 10.8 1.10 =16.6> 9120S 
E0-E3 12.0 1.13 =16.7 SO2ReSsS 
Bright Sa, Sap, Sab «611.9 0.90 —17.1 0.82 51 
SBO 11:4 1.10 —17.27 JO%53meeE 
SBb 11.2 1.10 -—17.4 0.86 26 
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tained empirically. Therefore the characterization of 
luminosity functions of particular types of galaxies is 
based on the mean M, of the absolute magnitude and on 
the dispersion oy. Table I gives the values obtained. 
The symbols a and @ in the first two columns of the 
table characterize the selection probabilities ¢:(m) of 
the various types. These were approximated by 


i) 


1 
$:(m) =—— (38) 


(On)! exp(— x?/2)dx. 


(m—a@)/B 


Thus, a stands for that apparent magnitude of galaxies 
of a specified type for which the probability of being 
included in the catalogue is equal to one-half. It is seen 
that the luminosity functions of several types of 
galaxies, some of them quite different morphologically, 
are rather similar. Because of this circumstance these 
particular types can be combined into groups, even 
though the relevant selection probabilities may be 
different. Thus the estimation of o? was performed 
separately for three groups of types of galaxies labeled 
faint, medium, and bright. The results are given in 
Table II. The last column gives the estimated standard 
errors of the estimates obtained. The last line gives the 
combined estimate (36) and its estimated standard error. 

Only the result for medium-bright galaxies is positive. 
The value of 14 681 obtained indicates the conditional 
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dispersion of redshift, supposed the same at all dis- 
tances equal to 121 km/sec. The other three estimates 
are negative. 

Even if one gives complete credence to the classi- 
fication of galaxies and to the groupings we used, these 
negative signs should not be surprising because of the 
large standard errors of the estimates. Taking into 
account the unpublished result of A. K. Kruszewski 
indicating that the error variance in measuring redshift 
is of the order of 2500= (50 km/sec)’, all the estimates 
obtained are compatible with the assumption that the 
intrinsic dispersion of redshift is negligible. The corre- 
sponding 95% confidence limit is approximately equal 
to 150 km/sec. Thus, the original estimate of Hubble 
appears somewhat exaggerated. 

All the above is written from the standpoint of the 
hypothesis of no error in our groupings of morphological 
types. However, we are aware that our groupings are 
likely to be faulty. In addition, as will be explained 
below, there is the possibility that, as a result of still 
another factor, the expansion of the universe, the 
figures given in Table II underestimate the true 
value of o°. 

Figures 1, 2 and 3 illustrate the fit of the regression 
lines to data for faint, medium, and bright galaxies. 
Both the empirical values of Y and the fitted expecta- 
tions of Y, given m, are plotted not against «?=10?/* 
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but against m itself. As a result, the regression lines are 
represented by exponential curves, with the coefficient 
of the exponential and the intercept determined by 
least squares as indicated above. 

It will be seen that the difference between the three 
regressions, ascribable to differences in the luminosity 
functions of the three categories of galaxies, appear so 
small that it is hard to see them, except in Fig. 4 where 
they are redrawn together to emphasize the distinction. 
However, this apparent similarity of regressions is 
deceptive. It is caused in the graph by a very large 
contraction of the vertical scale enforced by dealing 
with the squares of numbers representing velocities in 
thousands of kilometers per second. As a result, there 
are considerable concentrations of observational points 
near the horizontal axis and, in order to make these 
points visible at all, some of them are drawn im- 
mediately below the horizontal axis. Actually, they 
should be above this axis, in its immediate vicinity. 


Taste II. Estimates of variance of redshift of field galaxies. 


Group of types o2 (km/sec)? St (km/sec)? 


Faint —31 613=— (178)? 29 884= (173)? 
Medium +14681= (121)? 17 683= (133)? 
Bright —81 643= — (286)? 50 336 = (224)? 
Combined estimate — 4386 = — (66)? 14 567= (121)? 
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Fic. 3. Scatter diagram and | 
fitted regression line of Y on m. | 
Bright galaxies. i 


For small values of m there are no observational 
points and the three regression lines are so close together 
that they are indistinguishable. 

With these details in mind it is now easy to see that 
an accurate estimation of the variance o? of redshift is a 
problem of considerable delicacy. Small as the distances 
between the three regression lines are, it is obvious that 
an inappropriate classification of galaxies on which the 
regressions are based may shift the fitted curves either up. 
or down and, small as this shift may be, it may seriously 
affect the value of the estimate. 

Because of this fact, the appearance of negative 
estimates of an essentially nonnegative quantity should 
not be surprising. The system of types used by Huma- 
son, Mayall, and Sandage, essentially that of Hubble, is 
now being revised. In addition, of course, the combina- 
tion of types used in our calculations may easily be 
faulty. It may be hoped that, with the increase of 
observational material and a refinement of classification 
the method deduced may bring more reliable esti- 
mates of o?. 

Before concluding, we wish to point out that the 
negative estimates of o? may be an indication that our 
assumption of stationarity in time of the stochastic 
process representing the universe is incorrect. Although 
no numerical calculations are completed as yet, the 
formulae deduced indicate that the expansion of the 
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I’'1c. 4. Comparison of fitted regression lines of Y on m 
for faint, medium, and bright galaxies. 


universe combined with the absence of “‘steady state” 
implies that the regression line of Y on 2? is not linear 

but is convex upward. Asa result, if a straight regression 

of Y and 2? is forced on the data, then the fitted regres- 
sion line will have its intercept somewhat less than its 

true value and, when this is small, the estimate may 
_ easily tend to be negative. The details of this theory will 
_ be considered in another paper. 

Postscript. We are grateful to Professor Thornton L. 
Page for pointing out a reason why the assumption 
that the variance o? is independent of the morphological 
type of galaxies may be unrealistic. According to Page, 
the dispersion in redshift from Hubble law must be 
consistent with the equipartition of energy among 
galaxies. Since according to his findings (Page 1961a) 

_ the ellipticals are about 30 times as massive as the 
spirals, the latter may be expected to have a corre- 
spondingly higher o. This is an interesting possibility and 
is certainly worth investigating. Should it happen that 
does depend upon the morphological type, then formula 
(36) will yield an unbiased estimate of an average value 
of the corresponding variances. 
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The study of the dependence of o on the morpho- 
logical types is combined with practical difficulties. As 
indicated in Table II, the standard errors of the esti- 
mates found are very large, all of them of the same 
order of magnitude as the corresponding estimate. Thus, 
in order to be able to detect differences connected with 
morphological types, something must be done to 
improve the precision of the estimates. The requisite 
increase in the number of observations is probably 
prohibitive. However, it is just possible that the 
groupings of morphological types are rather seriously 
at fault. If so, then the improvement in this respect will 
decrease the residual variance in the regression analysis 
performed sufficiently to bring out the differences in o? 
if they exist. 
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NOTICE 
Request for Errata in Star Catalogues of Position 


HE United States Naval Observatory and the Yale University Observatory have been collaborating on the 
preparation of IBM punched cards of Zone Catalogues and certain other widely used catalogues. In the 
course of this undertaking numerous errors in the published catalogues have been discovered. A compilation of 
corrections is to be published soon after the completion of the card catalogues. 
The following are among the catalogues completed or in preparation: 
Yale Zone Catalogues, Trans. Yale Univ. Obs., Vols. 11-14, 16-27. 
Catalogue of Bright Stars (Schlesinger and Jenkins, 1940). 
General Catalogue of Trigonometric Stellar Parallaxes (Jenkins, 1952). 
AGK2, Vols. 1-15. 
Cape Photographic Catalogues 
Zones —30° to —40°, Ann. Cape Obs., Vols. 17, 18. on 
Zone —40° to —52°, Cape Zone Catalogue of 20843 stars, equinox 1900 (London, 1923). Proper motions of the 
stars in the Zone Catalogue of 1900. 
Zones —52° to —64°, Ann. Cape Obs., Vols. 19-20. 
N30, Astron. Papers Am. Ephemeris and Nautical Almanac, Vol. 13, Part 3, 1952. 
In order that the compilation of corrections may be as complete as possible, the authors request that other user: 


of these catalogues report any errors they have found to the undersigned. 
Dorrit HOFFLEIT 
Vale University Observatory 
Box 2023, Vale Station 
New Haven, Connecticut 
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NEW BOOKS RECEIVED 


Annuaire pour l’An 1961 (with Supplement for 1962), published by the Bureau of Longitudes, Paris. Price, paper 
back 40 NF, Composition 45 NF. 

Although containing in abbreviated form the principal astronomical tables of the standard national ephemerides 
this useful work has, in addition, an extraordinary variety of miscellaneous information, including, for example 
calendars ranging from the French Revolutionary to the Saints Days for 1961, and from Coptic to Cambodian 
Geodesy, refraction, tides, terrestrial magnetism, minor planets, vital statistics—all these and very much more ar‘ 
to be found in notes, tables, and graphs. Excellent review articles on such subjects as meteors, comets, stella 
parallax, novae and supernovae, artificial satellites, and a large French telescope make up more than half the book 
Use for reference, fun for browsing, very well indexed. (In French.) 


Aslronomischer Jahresbericht, 58. Band, published by the Astronomisches Rechen-Institut in Heidelberg. Price 
60-DM. 
Fifty-eighth annual volume of the unique and indispensable reference summary, covering in highly organizec 
and indexed form all astronomical-literature published in 1958. 


The Scientific American Book of Projects for the Amateur Scientist, by C. L. Strong. Published by Simon and Schuster 
New York. Price $5.95. 
The first section, on astronomy, contains nine articles of interest to teachers and astronomical gadgeteers 
Observational astronomers may find useful the detailed description of a transistorized drive and of the Leighto1 
electronic seeing compensator. 


Aussergalaktische Sternsysteme und Struktur der Welt im Grossen, by Henrich Vogt. Published by the Akademischi 
Verlagsgesellschaft, Leipzig, 1960. Price DM-20. (Volume 32 of Probleme der Kosmischen Physik.) 
Short review monograph on galaxies. A 70-page summary of some of the better-known observational materia 

is followed by succinct semiquantitative accounts of Newtonian, relativistic, kinematical, steady-state, ani 

variation-of-constants cosmologies. (In German.) 
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